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We present a comprehensive formalism for the description of primordial black hole formation in 
spherical symmetry based on the formalisms of Misner, Sharp, and Hernandez, which can be used 
to predict whether or not a black hole will form, and extract the resulting black hole mass when 
formation does occur. Rigorous derivations of all aspects of the formalism are provided, including a 
thorough investigation of appropriate initial and boundary conditions. We connect our formalism 


with numerous other approaches in the literature, 
are provided. We include animations of simulated 
material. 
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I. INTRODUCTION 

The very early universe is understood to consist of a 
generally isotropic and homogeneous plasma with small 
perturbations. These perturbations, thought to originate 
in quantum fluctuations driven outside the horizon during 
inflation mm, form the seeds of structure that give rise 
to the cosmic microwave background and ultimately the 
universe as we know it today. 

If the primordial perturbations are sufficiently large, 
it is possible that they might collapse under their own 
gravity, forming primordial black holes Eng. This re¬ 
quires perturbations to grow to nonlinear scales. The 
standard inflationary tale predicts Gaussian scale invari¬ 
ant perturbations with amplitudes ~ 10“®, meaning that 
the probability of forming a nonlinear perturbation is neg¬ 
ligible; indeed, primordial black hole production rates are 
predicted to be vanishingly small in this picture [S]. It is 
conceivable, however, that there exist features in the infla¬ 
tionary power spectrum at small wavelengths, well below 
cosmic microwave background radiation scales, that may 
allow perturbations to become sufficiently nonlinear as to 
grow under their own gravity and undergo gravitational 
collapse [SHS] (alternatively, large non-Gaussianities or 
isocurvature may also provide the necessary conditions). 

The idea of primordial black holes, whilst currently 
strictly hypothetical, is intriguing for a number of reasons. 
Astrophysically, an appropriate spectrum of primordial 
black holes may provide the seeds for the supermassive 
black holes found in the centers of galaxies [S] . An entirely 
different spectrum consisting primarily of lunar mass 
black holes may be able to explain cold dark matter cni- 
Smaller black holes are also of theoretical interest as a 
window into the process of black hole evaporation through 
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Hawking radiation m- Finally, should primordial black 
holes be identified, this would give insight into the small- 
scale structure of the very early universe, well beyond 
what the GMB and large scale structure can provide. 

The concept of primordial black holes was first intro¬ 
duced by Zel’dovich and Novikov [T^] in 1966, and again 
by Hawking in 1971 m- Theoretical work in the 1970’s 
identified rough criteria for primordial black holes to form 
mun], and early numerical hydrodynamic simulations 
were undertaken later that decade [TMTB] . Later work 
investigated the possibility of primordial black hole for¬ 
mation through exotic processes such as bubble collisions 
m, cosmic strings [201EU and phase transitions [22j . 

A large amount of work has gone into computing the 
number density and mass spectrum of primordial black 
holes from inflationary models laiTlIHlEs], and there 
exists an extensive literature constraining the mass spec¬ 
trum through astronomical, astrophysical and cosmolog¬ 
ical measurements [2M2^ (see also [301132] for reviews). 
The absence of primordial black holes has in turn been 
used to place constraints on other cosmological processes. 

In this paper, we are concerned with the formation of 
primordial black holes through inflationary perturbations 
in the primordial plasma. Most work in the literature 
deals with spherically symmetric perturbations, and we 
follow this trend. This simplification is justified by the 
result from peaks theory [33] that the largest peaks re¬ 
sulting from an appropriate probability distribution have 
a very strong tendency to be (almost) spherical. 

In order to understand the primordial black hole num¬ 
ber density and mass spectrum, one needs to know the 
threshold under which perturbations will form primordial 
black holes. Early work by Carr [T4l[Il] considered a col¬ 
lapsing region described by a closed Friedmann-Robertson- 
Walker (FRW) spacetime surrounded by a flat FRW uni¬ 
verse. In the radiation-dominated era, this was found to 
lead to a threshold value for the perturbation amplitude 
6c ~ 1/3, where 6 is defined to be the fractional mass 
excess inside the cosmological horizon when the overdense 
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region enters the horizon. 

Early hydrodynamical simulations of primordial black 
hole formation muH] were based on the Misner-Sharp 
formalism [33], which describes the gravitational collapse 
of a thermodynamic fluid under spherical symmetry. This 
formalism cannot evolve past the formation of a black 
hole, however. To our knowledge, the first numerical 
evolution capable of extracting the final state black hole 
mass was Niemayer and Jedamzik |35j . who modified stel¬ 
lar collapse code [33] based upon the Hernandez-Misner 
formalism which uses a null slicing condition to avoid 
singularity formation. Alternative formalisms have also 
been proposed [381141)] . A number of works have used 
these formalisms to investigate critical phenomena and 
shock formation near the threshold of black hole formation 


hole formation condition and appropriate cosmological 
boundary conditions. Section ^ addresses the subject of 
initial conditions in detail. We then derive the Hernandez- 
Misner formalism in Section jVj and cast it into variables 
more suited for cosmological evolution in Section VI We 
discuss various coding issues in Section jvn] Those who 
are only interested in seeing the supplemental animations 
may like to skip to Section |VIII[ where we demonstrate 
and explain simulations in our formalism. Finally, we 
conclude in Section |IX| A number of technical appendices 
are included. 


II. MISNER-SHARP FORMALISM 



Niemayer and Jedamzik found that the critical mass 
overdensity Sc ~ 0.7, much larger than had been previ¬ 
ously thought. However, Sasaki and Shibata (33] pointed 
out that this was due to their initial conditions containing 
an unphysical decaying mode. Musco et al. |45| investi¬ 
gated this further, arriving at a more modest 6c ~ 0.4. 
This issue gave rise to investigations of how to construct 
appropriate initial conditions for numerical simulations. A 
particularly nice method has been presented by Polnarev 
et al. [431 ITT] . 

While the formation condition for a primordial black 
hole is known to be roughly Sc 0.4, the precise value is 
dependent upon the density profile of the perturbation 
[48] . Polnarev et al. [49] [50] have recently constructed a 
new criteria which aims to capture the effect of the per¬ 
turbation profile dependence on the black hole formation 
condition. 

The purpose of this paper is to present a comprehensive 
formalism for the numerical evolution of spherically sym¬ 
metric perturbations in the early universe under the influ¬ 
ence of their own gravity. We employ the Misner-Sharp 
[84] and Hernandez-Misner m formalisms, building upon 
the ideas of Niemayer and Jedamzik |35j and Polnarev et 
al. [liSZ]. Our primary goal is to establish a clean and 
precise formalism from a somewhat murky literature. 

Our numerical formalism, while building on previous 
results, is entirely new. We correct errors in the litera¬ 
ture and present new techniques to improve accuracy in 
simulations. We analyze the regime of validity of various 
approximations, and provide a concrete description for 
how to connect inflation to initial conditions for numerical 
evolution. We have taken pains to be as rigorous and 
complete as possible, presenting our derivations in a ped¬ 
agogical manner. We do not give details of a particular 
numerical implementation, but point to issues that can 
lead to numerical troubles where appropriate. 

This paper is structured as follows. We begin by de¬ 
riving the Misner-Sharp formalism in detail in Section 
[Hi In Section [ml we apply the Misner-Sharp formalism 
to cosmology, and find variables more appropriate for 
this evolution. We pay particular attention to the black 


We begin by investigating Einstein’s equations for a 
spacetime containing a perfect fluid under the assumption 
of spherical symmetry. We follow the formalism first 
laid out by Misner and Sharp [34], and aim to present a 
detailed yet succinct derivation of the equations of motion. 

A spacetime is spherically symmetric if it possesses an 
50(3) isometry (3D rotation) group, where the action of 
this group on any given point (the ‘orbits’ of the group) 
are two-dimensional spheres m- The spacetime metric 
induces a two-dimensional metric on each such sphere, 
which must be proportional to the metric of a unit two- 
sphere, 

= d0^ + sin^ 6d(jP (I) 

for two coordinates 6 and (/>. This two-dimensional metric 
possesses three Killing vectors, associated with rotations 
around three axes. Demanding that the spacetime metric 
also possess these Killing vectors constrains it to take the 
form 

ds^ = gABdx^dx^ + R^{x^)dn^ ( 2 ) 

where x^ and x^ are arbitrary coordinates and the two- 
dimensional metric gAB depends only on and x^. 

The remaining gauge freedom in the metric ([^ lies in 
redefinitions of x^ and x^ . Common gauge choices are to 
choose coordinates such that = (a;^)^, known as the 
radial gauge, or to fix {x^)‘^gii = R^, known as isotropic 
radial coordinates. The gauge choice we employ is to fix 
goi — 0, which makes the metric diagonal. We choose 
coordinates t and A, and write the metric as 

ds^ = + e^dA^ + R^dfl^ . (3) 

Here (j)^ A, and R are functions of A and t alone, and we 
take i? to be a monotonically increasing function of A. 
The remaining gauge freedom lies in transformations of 
the form t = t{i) and A = A{A). 

We describe the matter content of the spacetime as a 
perfect fluid with stress-energy tensor 

= {p + P)uy,Uu + Pg^j,^ (4) 
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where p is the energy density, P is the pressure, and is 
the four-velocity of the fluid. Note that because this tensor 
is diagonal in the local rest frame of the fluid, it cannot 
describe the energy flow associated with radiation [M] , 
In particular, this means that free-streaming neutrinos 
cannot be described in this formalism. We discuss the use 
of a perfect-fluid description of the early universe photon 
bath in Appendix]^ It is convenient to use coordinates 
such that the coordinates move with the fluid, known 
as comoving or Lagrangian coordinates. The fluid four- 
velocity is then given by 

u* = e~^, M* = 0; i = r,9,4'. (5) 

The preliminaries complete, we now compute the equa¬ 
tions of motion. We begin with the Einstein equation 
where we work with G = c = 1. The 
time-time, radius-radius, and time-radius components of 
the equation are as follows, where we use primes and dots 
to denote derivatives with respect to A and t respectively. 


^ {e» + 


if-e 






R 


(e^RX + e^'^’R'X' - 2e^^R"^ 


Sire^P = \(R'f + 

R^ L 


0 = 


- + if - 2RR<j) + 2RR) 

XR' + 2R(j)' - 2R' 


R 


(6a) 


(6b) 


(6c) 


Except for the 69 and (j)(j) components, all other compo¬ 
nents of the Einstein equation vanish identically. The 99 
and (j)(j) components are related by symmetry, and as the 
Einstein equation and the conservation of stress-energy 
equation are related by the Bianchi identity, these angu¬ 
lar components can be neglected in favor of the simpler 
equations arising from conservation of stress-energy. 

We now use conservation of stress-energy, = 0. 

The t and r components are 


velocity, and therefore can exceed unity. The quantity E is 
simply an alternative variable to A, and must be positive 
due to the monotonic nature of R. The quantity m is the 
Misner-Sharp mass, a measure of the gravitational mass 
contained within the radius A. It must also be monotonic 
in A. The following related expressions are useful. 


e-^R = U + U<j) 

(9a) 

A/2^v ^ p' 1 py 

(9b) 

m' = AirPf R! p 

(9c) 


We now write the Hve equations of motion using these 
definitions. Begin by solving Eq. (7b I for </>'. 


P' 

p + P 


Next, write Eq. (6c) as 


XR' = 2{R! - R(f )'), 
2e'^U' 


X = 


R' 


( 10 ) 

(11a) 

(lib) 


Take Eq. (6a) and substitute in Eqs. (8a), (8b), (9b) and 
(11b) to obtain 

SirpR^R' = ^ ((1 -h C/2 _ r2)i?) . (12) 

Integrating this from 0 to A, imposing the boundary 
condition R{t,0) = 0, and using Eq. (8c) we find 


2 to 

r2 = 1 -h [/2-. 

R 


(13) 


Now , multiply Eq. ([7^ by 27tRR' and substitute Eq. 

to 

d 


(11aI. Some rearrangement leads to 
d 


0 = ^ {AnpR^R') + Att— (^PR^R^ . (14) 

Integrate this over the radial coordinate from 0 to A, then 
make use of Eqs. (8a) and (8c) to obtain 


Q^2Rp+{P + p){AR + XR) (7a) 

0 = P' -h (p + P)(t>' (7b) 

while the 9 and (/) components vanish identically. 

The following definitions allow us to bring the equations 
of motion into a more convenient form. 


U = e-‘^R 
r = e-^/'^R! 

pA 


m = Att 


pR^R'dA 


(8a) 

(8b) 

(8c) 


TO = -c^AttR^PU (15) 

where we once again employ the boundary condition 


R{t,0) = 0. The last equation is (6b). Substituting 
Eqs. (8a), (8b) and (9a) yields the fol 


lowing. 


2e-‘^RU = T^ -1-U^ + 2r2p^ - 8ttPR^ 


(16) 


This can be further simplified using Eqs. (10) and (13). 

P' m + A'kPR? a 


U = -e'f 


R'{p + P) 


+ 


P2 


V 


(17) 


U is the coordinate velocity of a fluid element, as the fluid 
element at comoving radius A is positioned at radius R. 
Note that this is only a coordinate velocity, not a physical 


We now have the six equations that constitute the 
Misner-Sharp formalism |34j . Three equations are evolu¬ 
tion equations for P, m and U (Eqs. (8a), (15) and 0), 
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and three equations are constraint equations for p, (j) and 
r (Eqs. ( [^ , (10), and (Hi)- Together, they express the 
content of Einstein’s equations in a convenient manner 
for numerical implementation. 


This system of equations must be supplemented by 
an equation of state P. In the original Misner-Sharp 
formalism, the equation of state was taken to depend on 
energy density and specific energy (kinetic energy density), 
and so thermodynamic relations were also required to close 
the system of equations. For early universe cosmology 
purposes, we can take pressure to be a function of ener^ 
density only, and so we only need the equation of stat^ 
Because it describes fluid flows, the Misner-Sharp for¬ 
malism has the capacity to develop shocks. A typical way 
to combat this is to introduce an artificial viscosity term 
into the equations in order to smooth out any shocks [52] . 
The effect of adding a a viscous term is to let pressure 
P = Po + Q, where Pq is the physical pressure, and Q is 
the artificial viscosity. For the moment we simply include 
it as an unknown artificial viscosity in the equations of 
motion, and leave a detailed discussion for Section [Vllj 
For the physical pressure, we will use Pq = wp, where w 
is the (possibly time-dependent) equation of state. For 
radiation, w = 1/3, as we detail in Appendix [A] 

We complete the formalism by specifying boundary 
conditions. We utilized R{t, 0) = 0 in the derivation of 
the equations. By the definitions of U and m, we then 
have U{t,0) = m(t, 0) = 0. By spherical symmetry, we 
have p'(t, 0) = P'{t,0) = 0 and thus (/'(t, 0) = 0 also. 
Using L’Hopital’s Rule we find that 


nience. 


Ue'^ 

(I9a) 

-P^A-kR^PU 

(I9b) 

e^f +^+4.Rp) 

(r'(p + P)^ R^^ J 

(19c) 

m! 

AttR^R' 

(I9d) 

P' 

p + P 

(I9e) 

Jt 

(I9f) 


The boundary conditions are R{t,0) = 0, m(t, 0) = 0, 
U{t,0) — 0 and r(t, 0) = 1, with an arbitrary (nonzero) 
boundary condition prefit). For later convenience, we 
also include the following expression. 

p = + + ( 20 ) 


When artificial viscosity is turned off and w is constant, 
we can integrate Eq. (19e) to obtain 


= 


Pref 


w / (1+it;) 


( 21 ) 


where pref is some reference density which is simply re¬ 
lated to the boundary condition on (j). This formula is 
useful for analytic analysis. 


lim = 1 
A-J-O 


lim 

A -)-0 


SirpR^R' 

R' 


= 1 . 


(18) III. COSMOLOGICAL APPLICATION OF THE 
MISNER-SHARP FORMALISM 


The only remaining boundary condition is to specify (j) at 
one of the boundaries. This is simply a gauge condition 
on the time coordinate, and is usually set by matching to 
asymptotic behavior. 

Initial conditions must be set for each of the evolution 
variables, namely m{to,A), U{to,A) and R{to,A). The 
first two of these correspond to the local energy density 
and fluid velocity, while the third of these is the initial 
gauge condition. The fluid provides the only dynamical 
degree of freedom; the gauge condition simply evolves 
with time as dictated by the fluid. 

We summarize the Misner-Sharp equations for conve- 


The Misner-Sharp formalism was originally designed 
to study stellar collapse, where pressure vanishes after 
some radius, and the metric smoothly connects to the 
Schwarzschild metric. In cosmology, neither pressure nor 
energy density vanish at large radii, and spacetime is 
not asymptotically flat. Thus, particular care must be 
given to the boundary conditions in the Misner-Sharp 
formalism in this case. In this section, we use the Misner- 
Sharp formalism to describe cosmological evolution. We 
begin by investigating the Friedmann-Robertson-Walker 
(FRW) limit of the equations, before looking beyond the 
background evolution to find variables more suited to 
cosmological evolution. 


A. Background Evolution 


A number of papers in the primordial black hole literature employ 
thermodynamic relations in their system of equations despite this. 
To our knowledge, this is unnecessary. 


We begin by exploring the background cosmological 
evolution. This requires that the energy density is con¬ 
stant in space, yielding p/ = 0, where the subscript b 
indicates the background FRW value. This in turn leads 
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to = 0 from Eq. (19e) where we ignore artificial viscos¬ 
ity. By choosing our boundary condition appropriately, 
we c an set (pb = 0. This gives us p^ef = Pb- Equation 
(19d) can now be integrated to obtain 


The metric then becomes the curved FRW metric, 


ds^ = —dt^ + a{ty 


dA^ 

l-kA^ 


■ A^dfl- 


(34) 


47r 3 

mb = -^PbRb ■ 


( 22 ) 


Using these results, the evolution equations (19a), (19b) 


and (19c) become the following. 


Rb = Ub (23) 

Pb =-3^Pb{l + w). (24) 

rib 

• 47r 

Ub = — —pbRb{^ + 3w). (25) 


where we used Thus, the Misner-Sharp 

equations yield the FRW cosmological evolution, as ex¬ 
pected. 

Let us specialize to the case of A: = 0 (flat FRW) with 
constant equation of state w. As we aim to model cosmo¬ 
logical evolution in the radiation dominated era immedi¬ 
ately after inflation has suppressed all spatial curvature, 
this is a reasonable regime to investigate. 

We can then solve the background cosmological evolu¬ 
tion as follows. The continuity equation is integrated to 
obtain 


In particular, note that there are no spatial derivatives 
here, and so the evolution of each function will be in¬ 
dependent of A. In particular, this allows us to write 
Rb{A,t) = a{t)A where we exploit reparameterization 
invariance to write Rb{A,to) — A for our initial data, 
selecting a{to) = 1. The function a{t) is the usual scale 
factor of FRW cosmology. Substituting this description of 
R into the equations of motion then yields the following. 

Pb =-3-pbil + w) (26) 

a 

Ub = Ad (27) 

• Air 

Ub =-—pbaA{l+ 3w) (28) 

The first of these is recognizable as the continuity equation 
in cosmology. The last two of these combine to yield 

- = -^Pb{l + 3w) (29) 

d O 


Pb = poo- 


3(l-|-iu) 


(35) 


where po is the FRW energy density at time to- The 
Friedmann equation then becomes 


d = 



(36) 


where we have defined a = 2/3(l-|-w) for convenience 
{a = 1/2 for w = 1/3). This can be integrated to obtain 


a(a'/“-l) = /^(t-Ao) (37) 

where the boundary condition a{to) = 1 has been used, 
and we assume w y —1. It is convenient to set to by 


Sttpo 

such that 


(38) 


which is the acceleration equation. Combined with the 
continuity equation, this can be used to obtain 



(39) 


9t(d^) = ^dtipbo^) 


which integrates to give 
2 


Stt 


= H^ = ^pb - 


(30) 


The Hubble parameter is then given by 


H = 


a 

t 


(40) 


(31) 


B. Beyond the Background 


for some constant of integration k. This is the Friedmann 
equation with spatial curvature k. We can multiply this 
by Rl and substitute for Ub to obtain 

= kA‘ (32) 

Kb 

and thus 

Tl = l-kA‘^. (33) 


While the full Misner-Sharp equations describe grav¬ 
itational collapse in cosmology, they are not written in 
terms of the most efficient quantities, particularly from a 
numerical perspective. As we can solve the background 
evolution exactly, it is useful to factor out this aspect 
of the solution from the variables, in order to enhance 
accuracy and stability. In this subsection, we transform 
the Misner-Sharp equations in such a manner. This sub¬ 
section largely draws upon and improves the ideas of 
Polnarev, Nakama and Yokoyama m- 
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R = aAR — RbR 

(41a) 

P = PbP 

(41b) 

47r n _ 


m = —pbR fh 

3 

(41c) 

U = HRU 

(41d) 

P = PbP = pbp{w + Q) 

(41e) 


where the tilded quantities are unity for FRW evolution, 
and are even functions (Q, which represents a dimen¬ 
sionless form of the artificial viscosity, is also an even 
function but vanishes for FRW evolution). The choice to 
scale fa and ?7 by i? instead of Rb is made because the 
expected local quantities depend on physical radius rather 
than the expected FRW radius. In more pragmatic terms, 
this choice of quantity simplihes the expressions that fol¬ 
low. The Misner-Sharp equations of motion become the 
following in this notation. 


P' 

p + P 


(Ue^ - l) 


-R = R 

. , AR _, 

p = m -\ -=—m 

3{ARy 

—m = —m — 3Ue^ [P -\-rh] 
Ha V / 

= 1 -H H^R^ {u^ - m'j 


iJ a 


F^ 


P' 


H'^RR'{p + P) 

+ i {2U‘^ +fh + W 


(42a) 

(42b) 

(42c) 

(42d) 

(42e) 

(42f) 


Here, we’ve made liberal use of the background expres¬ 
sions, the continuity equation and the Friedmann equa¬ 
tion. 


For the purposes of numerical evolution, the evolution 
slows down significantly as time increases due to the 1/t 
that appears in the equations of motion (as part of the 
Hubble parameter). To address this, it is useful to work 
in logarithmic time. Define 


^ = ln 



— Ina 
a 


(43) 


such that 

dt = -d^. (44) 

a 

We can express the Hubble and scale factors in terms of 
^ as follows. 

a = e“« (45) 

H = —(46) 
to Rh 


Here, we define Rh = to/o- as the horizon radius at the 
beginning of the evolution. For later use, it is convenient 
to note that 


HaRu = , 

the horizon scale evolves as 

Rh{P) = = RhS-^ , 

and that the Friedmann equation can be written as 


Stt 


= -2? 


^ “ 3 “ i?2 


(47) 

(48) 

(49) 


H 


It is also convenient to define A = ARh with A dimen¬ 
sionless. Using these variables and the definition (38) for 
to, the equations of motion become the following. 


P' 


d^R = aR {Ue’^ 


- 1 


- . AR ^, 

p = m-\ - —m 

Z{ARy 

d^fh = 2m — 3aUe^ (^P + fh^ 

r2 


F^ = 




H 


(50a) 
(50b) 
(50c) 
(50d) 

l-a)C ^2 


= 2 ( 1 - 


d^u = u-. 




P' 


AR{R + AR'){p + P) 


-f - (21/2 +fh + ?,P 


(50f) 


Here, we repurpose primes to denote derivatives with 
respect to A, and found it useful to use f in favor of F. 
When the artificial viscosity is vanishing, we have 

= p-3o.w/2 ( 5 ^) 


These equations are supplemented by the boundary condi¬ 
tions p' = 0, R' = 0, m' = 0 and U' = 0 at H = 0 (i.e., all 
functions are even functions of H). A boundary condition 
on (j) is also required; we suggest using the result ( |5T| ). 
Initial data consists of specifying R, U and fh (or p) at 
time ^ = 0. 

The values of m, R and U at the origin are only im¬ 
portant for computing derivatives, as the quantities m, 
R and U already vanish at the origin. If some evolution 
scheme that omits the origin while requiring fh, R and 
U to be even functions is used, then the origin can be 
ignored altogether. If not, then the functions should be 
evolved subject to the constraint that these functions have 
a vanishing derivative at the origin. We also point out 
that the evolution equation for U is not singular at the 
origin as P is an even function; with artificial viscosity 
vanishing, P = wp, and 


lim -j = p"(0) = -m"(0) 
A-j-o A 3 


( 52 ) 
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where we’ve used to' = 0 at the origin. 

The equations (501 are the complete dimensionless cos¬ 
mological Misner-Sharp equations of motion. The only 
time explicit time dependence is in the expression for T^, 
and there is no dependence upon po or to at all. This 
implies that for any initial to, R and U, the evolution is 
the same no matter the initial time or density. 


C. Black Hole Formation 


Thus, one should check for the conditions U < 0 and (57) 
as the conditions under which a black hole forms. 


After the apparent horizon forms, a singularity at the 
origin forms soon afterwards. The Misner-Sharp formal¬ 
ism cannot evolve past the formation of this singularity, 
and so a new technique is needed if more information 
about the resulting black hole (such as its mass) is de¬ 
sired. This brings us to the Hernandez-Misner formalism, 
which we introduce in Section Ivl 


We now possess the tools to evolve initial density per¬ 
turbations to investigate whether they form black holes or 
not. We thus need a method for detecting the appearance 
of a black hole. To do this, we search for the creation 
of an apparent horizon, which if detected, points to the 
existence of an event horizor0 We point the reader to §7 
of [S3] for the following results. 

Define to be the timelike unit vector orthogonal to 
slices of constant t, and s“ to be the spacelike unit vector 
orthogonal to slices of constant A. Together, these define 
an outgoing null vector = (n“-|-s“)/-\/2- In our metric, 
= (e-'^,e-^/2 0,0)/y2. 

The expansion of outgoing radial null geodesics is given 
by 0 = akb where to“^ = -|- rAn^ — s“s^ is the 

two-dimensional metric on the hypersurface of constant 
t and A. An apparent horizon forms when 0 < 0. It is 
known that 

0 = (53) 

for spherically symmetric spacetimes. For our metric, we 
have 


0 = ^(c/ + r). 


(54) 


As such, an apparent horizon forms when {7 + T < 0. As 
T > 0, this requires U to be negative. Using Eq. (13), 
the condition of black hole formation becomes 


2m 

li 


> I 


(55) 


which is recognizable as the amount of mass inside an 
circumferential radius R smaller than the Schwarzschild 
radiu s, and thus satisfying the hoop conjecture. Using 
(41c) and the Friedmann equation, this yields 


H^R^m > 1. 


Written in terms of our variables, this becomes 


(56) 


(57) 


^ Note that the lack of an apparent horizon does not necessarily 
imply that there is no event horizon, but this is typically limited 
to specific spacetime slicings. For our purposes, we may assume 
that an apparent horizon will form before the central singularity. 


D. Outer Boundary Conditions 


Let us now consider a situation where we wish to inves¬ 
tigate gravitational collapse near the origin, and desire to 
smoothly connect to an FRW evolution at some radius. 
Suppose that our region of interest is A < Aq on the 
initial time slice. By investigation of the Misner-Sharp 
equations, we see that it is insufficient for the local quan¬ 
tities {R, U, p and (f)) to be equal to their FRW values to 
ensure FRW evolution for A > Ag, as an incorrect mass 
inside the radius Ag will then cause the evolution to differ 
from FRW. Thus, in order to ensure FRW evolution for 
A > Ag, we must also require TO(Ag) = TOb(Ag). This 
matching condition implies that any overdensity must be 
surrounded by a corresponding underdensity in order to 
smoothly connect to FRW. It is important to note that 
this condition is not guaranteed to be preserved under 
time evolution, as we discuss below. 

In terms of the tilde variables, the matching condition 
is that TO, U and R become unity at the boundary, and 
remain unity beyond it. If m{Ao) = TOb(Ag), then 


rAo 

0 = Att R^R'{p-pb)dA'. 
^0 


Letting p = 1 + <5^, this is equivalent to 


(58) 


0 = / R^R'pbSpdA'. 

Jo 


Let TO = 1 + 5m- Then from Eq. (50), we have 
1 


<5p = 


3(AR)'^{ARy L 


(Ai?)^(5„ 


(59) 


(60) 


Inserting this into Eq. ( |5^ along with R = aAHuR 
yields 


rAo r- 


0 = 


{ARf5m dA' = {ARf{rh - 1) 


A—AqI Rh 

(61) 


If instead of matching to FRW at Ag, we wish to match 
to FRW asymptotically, then our condition becomes 

Jim (Ai?)^(TO — 1) = 0 . (62) 

A—>-oo 















This extends the linearized argument given by Polnarev 
and Musco [46]. 

For the purposes of matching to FRW numerically, it is 
helpful to have the initial data approach the background 
values at the boundary of the computational domain. 
However, as the evolution equations depend on derivatives 
of the local variables, data that matches to FRW at the 
boundary at one time is not guaranteed to match at a 
later time. In particular, a central perturbed region tends 
to gradually expand out, particularly when it doesn’t 
form a black hole. 


This means that there is no physical boundary condi¬ 
tion at the outer boundary of the computational domain. 
However, in order to solve the equations, an outer bound¬ 
ary needs to be specified, as there are two characteristics 
to the equations of motion, one traveling inwards and one 
traveling outwards. Without a boundary condition, no 
control over the ingoing characteristic is exercised, and 
the evolution rapidly becomes unstable. Just like a one¬ 
dimensional wave equation on a finite domain, a single 
boundary condition is required at the outer boundary. 

A number of works in the literature have employed the 
Dirichlet boundary condition p{Aq) = pf,. This boundary 
condition does an excellent job of maintaining stability at 
the outer boundary, but is unphysical in that it asserts 
that what is outside the boundary remains FRW for all 
time, contrary to the argument given above. As such, it 
creates a reflecting boundary condition where outgoing 
waves reflect negatively off the outer boundary. We ex¬ 
perimented with a Neumann boundary condition p' = Q 
at the boundary, and as expected, found that waves were 
reflected positively. 

One way in which to ameliorate this issue is to begin 
the evolution with a sufficiently large domain such that 
traveling waves do not reach the boundary before it can 
be determined whether a black hole has formed/will not 
form. However, this is computationally expensive. One 
solution to this, advocated by Nakama et al. [49], is to 
begin with a reasonably small domain that matches to 
FRW on the boundary, and watch for the boundary to 
deviate from FRW by a small quantity. Then, expand 
the domain by appending FRW values to the edge of the 
domain, and watch the new boundary for deviations. This 
would seem to work well, except for the fact that they 
required pinning p at the boundary to its FRW value, and 
watching U deviate, whereas in principle, both should be 
allowed to vary. In a similar vein, Hawke and Stewart 
[32] suggested adding extra gridpointss at the boundary 
of the domain as the evolution proceeds. 


We therefore desire a boundary condition on the outer 
boundary that allows waves to travel past the outer bound¬ 
ary with minimal reflection. To do so, we turn to the 
method of characteristics, the detailed application of 
which can be found in Appendix]^ We extract inward (^ 2 ) 
and outward (ua) moving characteristics from solutions 
to th e linearized equations ( |70| ). The equations of motion 
(B9| that U 2 and satisfy are inhomogeneous, with the 


inward and outward moving characteristics sourcing each 
other. This complicates the problem of constructing a 
nonreflecting boundary condition. In particular, we can¬ 
not simply set the amplitude of inward waves to be zero at 
the boundary because they are sourced by outward waves 
that have passed the boundary. Thus, the condition on 
U 2 for a nonreflecting boundary depends on the entire 
history at the boundary. 


For reasons of practicality, we pursue a local condition 
which gives very small reflections. Following a suggestion 
by Kidder et al. [54] , we try a Robin boundary condition. 


(A"u 2 )^ = 0. This condition is an assumption that U 2 
behaves like ~ outsi de the computational domain. 
Combining this with Eq. (B9b I to compute a condition 
we are able to impose, we find 


d^U2+Cs^U2= ^Ui + ^{U2+U3), (63) 


where Cg = e^/^/-\/T2 is the linearized wave speed, and we 
have specialized to w = l/d Converting this back from 
the characteristic variables, we obtain 


d,u2=^-f-S'^-^ (64) 

4 cs 

and so the following boundary condition for 
1 2c^ 

d^6u = -^6m-^S'^ 

+ {d^Sm + Cs^m) + ~^^'m ~ (65) 

From experimentation, we find that this condition works 
well for n = 3. From further experimentation, we also 
found that n = 1 for the term and n = 3 for the 
CsS'j^ term worked even better. We thus suggest using 

d^Su = + - CsS'u 

( 66 ) 


as the outer boundary condition. 

In our experimentation, we found that the linear ap¬ 
proximation is typically a good approximation at the outer 
boundary. However, as we know the nonlinear speed of 
wave propagation for a perfect fluid (as is needed when 
computing the CFL condition for numerical evolution, see 
Section VH), we tried substituting this into our boundary 
condition instead of the linearized speed of sound. We 
found that this didn’t noticeably decrease the reflected 
portion of the wave. 

Having a minimally reflecting boundary helps with 
the numerical evolution in a number of ways. It allows 


® The speed of sound isn’t = 1/3 as one might expect. This is 
due to the use of dimensionless variables; converting back to t 
and A yields the usual result. 
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Figure 1. Demonstration of our boundary condition. This plots density as a function of comoving radius at various times. 
The top curve is the earliest curve in this set, and the bottom curve is the latest. The curve in purple has an outer boundary 
at A = 20, while the curve in teal has an outer boundary at ^ = 40. The y axis is logarithmic, but the overall scaling is 
unimportant; it was chosen this way simply to show data at different timesteps clearly. 

This run did not form a black hole, but instead sent a wave of material propagating outwards. We see that our boundary 
condition allows the pnrple curve to track the teal curve remarkably well, despite terminating halfway along the domain. The 
largest deviations occur in the bottom two curves, which show that there was a small reflected component at the boundary. As 
the characteristics evolve as ~ 1/^, even a very small reflected component grows as it moves towards the origin. 


waves to escape the computational domain, which in turn 
allows us to reduce our domain size without affecting the 
evolution. This significantly increases the computational 
speed of our evolutions. 

There are two caveats to this boundary condition. The 
first is that it was found by using linear perturbation 
theory. As such, if a wave with nonlinear amplitude hits 
the boundary, the condition doesn’t work nearly so well 
(although it works better than the reflecting boundary 
conditions). The second is that it is physically possible 
(though unlikely, for appropriate initial conditions) for a 
wave to exit the boundary and then fall back in again. Our 
boundary condition does not allow for this occurrence. 
Regardless, the boundary condition works remarkably 
well, as we show in Figure]^ 


IV. INITIAL CONDITIONS 

The formalism presented above has evolution equations 
for m, R and 17, and constraint equations for calculating 
p, r and (f). As such, initial data for m, R, and U should 
be specified on a surface of constant t (although m may 
be substituted for p or (p A desired). We note that the 
initial data for R is simply a gauge choice; the physical 
initial data is m{R) and U{R), describing the fluid energy 
density and velocity. 

It is intended that the initial conditions are initially 
generated from inflationary perturbations. Small per¬ 
turbations outside the horizon grow, enter the horizon 
and become nonlinear before undergoing gravitational col¬ 
lapse. While these perturbations are small, they are well 
described in the linear regime. Linear fluid perturbations 
in cosmology are known to have two modes: a growing 
mode, and a decaying mode. In the time between the end 
of inflation and the beginning of our numerical evolution. 
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the decaying mode should have had sufficient time to 
vanish almost completely. As such, the initial data should 
represent this, and be comprised of the growing mode 
on 130 Disregard for this issue affected early work on this 
subject 1^, as pointed out in Ref. [55] . 

More recently, [35] developed a formulation that selects 
initial data comprising only a growing mode in the linear 
regime. This formalism was expanded in m. and utilized 
in |49| . In this section, we clarify and further develop this 
formalism. 


A. Linearized Equations of Motion 


We begin by computing the linearized equations of 
motion. The quantities that we need to linearize are 
i?, TO, and U, which are all unity at zeroth order (for 
simplicity, we take the artificial viscosity Q = 0 at linear 
order; see below). For efficiency, we also linearize p and (j). 
We do not need to linearize T, as T^ only appears in the 
equations of motion as a coefficient of p', which is itself 
first order. 

We linearize these quantities as follows. 



R=1 + €6r 

(67a) 


U = l + e6u 

(67b) 


rh = 1 + e6m 

(67c) 


p = l + e6p 

(67d) 


e* = l + e6^ 

(67e) 

Here, we use e 

as an order counting parameter. 

It is 

straightforward 

tions. 

to linearize Eqs. (50) from these defini- 


j. 3wa 

<30 = ^ 6p 

( 68 a) 

= oi {6u + 6^) 

( 68 b) 

6p = 6m + ~^^'m 

( 68 c) 

d^6m = iaw6m — ‘26u — 26^ — 3aw6p 

( 68 d) 

O' 

d^6u = (1 - 2a)6u - <^0 - - ((5m + 5w6p) 


Swa^6'p 2 (l-a)f 

2 A 

( 68 e) 


Note that Sn only appears in the evolution equation for 
5r, and is otherwise absent from this system of equations. 
This justifies the particular choice of variables in Eqs. 


^ It is not critical that the initial data consist of only a growing 
mode if sufficient time is allowed for the decaying mode to decay 
before any measurements are taken from the simulation. However, 
it is convenient to remove the decaying mode in order to begin 
simulations closer to black hole formation, and also to better 
simulate initial conditions from inflation. 


(411. We can thus ignore 5r entirely for the moment, and 
focus on 6m and 6u- Eliminating (5^ and 5p, we obtain 
the following evolution equations. 


d^6m = Swadm — 26u (69a) 

d^6u = (1 — ‘2a)Su — —6m 

(69b) 

The first of these can be solved for 6u, which can then 
used to eliminate 6ij in the second equation. This leaves 
us with a single second-order PDE. 

d‘^6m - (3 - 5a)d^6m + [(1 - 2a)3tc - 1] a6m 

= (^^+C) (70) 


It turns out that this equation is separable, and can be 
solved using a mode expansion. Unfortunately, all of the 
spatial modes diverge as A —> 0 , indicating that the linear 
approximation breaks down for small A. We thus look 
towards other methods. 

We begin by neglecting the RHS, which reduces the 
equation to a straightforward ODE. The solution will be of 
the form <5^ = C'exp(/I^), which yields the characteristic 
equation 

/32-(3-5a)/3-b[(l-2a)3w-l]a = 0. (71) 

This has solutions 

w — 1 

/3i = ^ ^ ~ 

For — 1 < re < 1, the first solution is a decaying mode, 
while for w > —1/3, the second solution is a growing 
mode. For our purposes, we will assume w > —1/3, and 
take only the second solution. 

5m(Ae)= (73) 


Here, we let 6mo = <5m(A,^ = 0) = 6m{A,t = to). Using 
this in Eq. (69a), we can compute the solution for 6u. 


ct. 


Su{A,0 = -^SmiA,0 


(74) 


To compute 6a, we take Eq. ( 68 bI and substitute in the 
expressions for 6^, 6p, 6u and 6m- 

= oi{6u + 6(j,) (75) 

= - ^ [2(1 - a)6mQ + awA6'mo\ (76) 


The solution to this is 

= — -^ ^mO 


1 -I- 3w 


-A6L 


mO 


p2(l-a)4 


(77) 


where we have discarded the constant homogeneous solu¬ 
tion which is a part of the background evolution. 
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These expressions for 5m, 5u and Sr provide us with 
exactly what we were looking for: a relationship between 
the dynamical variables that picks out the growing mode 
in initial data. Note that the time dependence of all three 
quantities is the same. These results are equivalent to 
Eqs. (4.24), (4.25) and (4.28) in Ref. [47]. 


B. A Derivative Expansion 


We now come back to the terms we discarded on the 


RHS of Eq. (70). The solution we found for 5m in Eq. (73) 


is the zeroth order solution 5m in a derivative expansion 
solution to Eq. (701. Here, we flesh out the details of the 


derivative expansion. 

We would like terms with derivatives acting on them 
to be suppressed compared to terms without derivatives. 
In order to accomplish this, we formally let A Ajy/e, 
where e is an order counting parameter for our derivative 
expansion (as distinct from e, which is the order counting 
parameter for the perturbative expansion). Equation (70) 
becomes the following under this transformation. 


d^5m — (3 — ^a)d^5m + [(1 — 2a)3w — 1] a<5„ 
A 




(78) 


We see that this demotes the terms on the RHS to con¬ 
tribute at first order in the derivative expansion (our 
solution ( |7^ is first order in perturbations, but zeroth 
order in the derivative expansion). 


Let us write our solution to Eq. (78) as 




(79) 


n—0 


Let us then look at the time dependence of this equation 
of motion. At zeroth order, the solution has time depen¬ 
dence given by exp[2(l — a)^]. The first order equation 
of motion then becomes 


dpi - (3 - 5a)d^5l + [(1 - 2a)3w - 1] a5l 

= + • ( 81 ) 

The time dependence of the particular solution must then 
be given by exp[4(l — a)^]. The full solution is 


= e4(i-“)€ 


2w 


7 -f 30w -f 27^2 




Co 


(82) 


Assume that at n**' order, the time dependence of the 
solution is given by exp[2(n-|-l)(l —a)^]. Then at (n-|-l)‘*' 
order, the equation of motion is 


- (3 - 5a)d^5p^ + [(1 - 2a)5w - 1] 

' 4 ( 5 ”' 


n+1 




A 


-b(5” 


(83) 


where in the second line, we substituted the time depen¬ 
dence from the n**' order solution. The particular solution 
to this ODE will be proportional to exp[2(n-I- 2)(1 — a)^], 
and hence by induction, the time dependence at n**' order 
will be exp[2(n -|- 1)(1 — Q;)'C]- All spatial dependence at a 
given order is determined by the spatial dependence at 
the previous order alone. Thus, a recursion relation can 
be constructed relating the solution at n*** order to the 
solution at (n — 1)*^ order. 

We thus (in principle) have the full solution to Eq. (70) 
as an infinite sum in a derivative expansion. 


C. Beyond Linear Order 


where our order counting parameter keeps track of the 
order of each solution. Our desire is that terms of succes¬ 
sively higher order become less and less important. 


We can substitute this solution into Eq. (78), and 
demand that the equation of motion be satisfied order by 
order. At 0{A^) (n > 0), the equation of motion becomes 


dpi - (3 - 5a)d^5l + [(1 - 2a)Zw - 1] a5l 




A 


■ 5t-^^" 


(80) 


Having solved the equations of motion at first order in 
perturbations through the use of a derivative expansion, 
we now look to higher order in perturbation theory. For 
the moment, let us neglect all spatial derivatives in the 
equations of motion, and consider the effect of expanding 
the equations of motion to higher order in perturbations. 
Let us use the expansion 


X = l + J^e^X^(A,a (84) 

n—1 


The LHS here is identical to the LHS of Eq. (70), and so 


has a homogeneous solution identical to Eq. (73), which 
can be absorbed into 51- Only the particular solution is 
then of relevance. The RHS of this equation consists only 
of known functions from the previous order. As such, this 
is an ODE rather than a PDE: the spatial dependence 
is completely described by the spatial dependence of the 
previous order. 


where X represents rh, U, R, p and At 0(e"), the 
equations of motion will have homogeneous terms that 
depend on the functions A„, and source terms that depend 
on Xm with m < n. The homogeneous terms will have 
equations of motion identical to the first order equations 
(68), and thus, identical solutions, which can be absorbed 


into the first order solutions. We thus look at particular 
solutions only. 
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The first order solution has time dependence given by 
exp[2(l — a)^], and therefore, the source terms for the 
second order equations will have time dependence given 
by (exp[2(l — The solution to the second order 

equations will thus have time dependence exp[4(l — a)^]. 
Assuming that for 0(e") and lower, the solution have time 
dependence exp[2n(l — a)^], we can again use induction 
to show that the solutions at all orders have this time 
dependence. 

It is curious that the time-dependence of both the 
derivative expansion and the perturbative expansion in¬ 
crease with powers of exp[2(l — a)^]. This motivates 
combining the two expansions into a single expansion 
wherein each higher order has time dependence given by 
an increasing power of exp[2(l — a)^]. Indeed, this is 
achieved by the expansion 

OO 

X = l + J2 (85) 

n—1 


in order to obtain higher order solutions from lower order 
ones. We do not aim to rederive their recursion relations. 
Instead, we provide expressions to assist in the construc¬ 
tion of initial data, where we specialize to w = 1/3. 

We provide expressions for to, U, R, p and e'^ to second 
order in perturbations, based of! the function Smo- We 
expand these quantities as follows. 


0 = 1 + e^Sl^iA) + e^^Sl,{A) 

(86a) 

U (A, 0 = 1 + ^^SIj{A) + e^^5lj{A) 

(86b) 

R{A, 0 = 1 + e^<5)j(A) -1- e^^(5|j(A) 

(86c) 

P(A, 0 = 1 + e^Sl{A) + e^«(5^(A) 

(86d) 

(eO(A,0 = l + ef<5^(A) + e2«(5^(A) 

(86e) 

first order solutions are 



where X represents m, U, R, p and e'^, and we have let 
e = e. The zeroth order terms are unity, while the first 
order terms are first order in perturbations and zeroth 
order in the derivative expansion (there are no terms 
that are zeroth order in perturbations and first order in 
the derivative expansion, as unperturbed terms have no 
spatial dependence). Second order terms are the sum 
of second order terms in perturbations and terms that 
are first order in both perturbations and the derivative 
expansion. 


The equations of motion in this expansion are as in Eqs. 


(50) except for the d^U equation, which upon substituting 
and —>■ multiplies p’ by e. As 

previously, e is only an order counting parameter. 


<5^ = (87a) 

5u = -^5m0 (87b) 

= + (87c) 

<5p = Smo + (87d) 

(87e) 


D. Second Order Solutions 


The expansion (85) is essentially the same as was used 
in Ref. j47j . where they constructed recursion relations 


where i5mo(A) is a free function. These expressions are 
equivalent to Eqs. (73), ( |74| ), ( [77| , ( |68a[ ) and (68c). 

The second order solutions are as follows. 


e = f (K - 6^ - <5^) + ^ - 3J^) + ^ 


S'^ - — 
'^^"20 


<5i= TV 


Si, (SL + Sl- 2Sh) - ^ 


SI 

2A 


16 


‘^SrSI/ -b 4(5// — S^ + Sl I -Sl — Sh — Sji 


Si = Si + A 


A2' 

-f-isl-sM 


si = ^, mr - 851] 


(88a) 

(88b) 

(88c) 

(88d) 

(88e) 
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The contributions at second order in perturbations (terms 
containing i5^) and first order in derivative expansions 
(terms containing S') are clearly evident in S'^ and Sfj. 
These expressions allow for the construction of initial 
data to second order that consist solely of growing modes, 
based on the choice of a single function Smo)S\S). 


E. Constraints on Initial Data 


There are a number of constraints on initial data, arising 
from both physical and desirable considerations. Here, we 
construct the various constraints, and then apply them to 
the solutions found above. We assume that initial data is 
specified in terms of m, R and U. 


1. Physical Constraints 


The most basic physical constraint is that m and R are 
positive, as neither m nor R can become negative. As U 
can be negative, so may U. By spherical symmetry, we 
also demand that m, R and U are even functions, which 
specifies that they have zero derivative at the origin. Note 
that this precludes the possibility of a sharp cusp in the 
energy density. 

Next, we require both R and m to be monotonically 
increasing with A. These yield the constraints 


0<R + AR' 


0 < m + 


AR 

—= —m 

3{R + AR') 


(89) 

(90) 


respectively. The second of these is equivalent to p > 0. 


From Eq. (50e), we have 

R^ (m - C/2) 


< 


A2 


(91) 


by requiring f ^ > 0 on the initial slice. This is equivalent 
to requiring that the metric component gAA is positive, 
so that we have the correct metric signature. 


A similar requirement arises from demanding that the 
initial data does not already describe a black hole. From 
Eq. (57), this requires 


R^rh < 


1 

l2 


(92) 


whenever U is negative. 


8. Desirable Constraints 


In addition to the physical constraints described above, 
there are a few desirable constraints also. The first is a nu¬ 
merical issue: it is recommended that initial data connect 


smoothly to FRW at the boundary of the computational 
domain. This improves the accuracy of the evolution by 
reducing unphysical reflections from the boundary. 


The second desirable constraint, as discussed previously, 
is that the initial data should consist of a growing mode 
only. This motivates us to use the solutions found in the 
previous subsection, which have been crafted to consist 
purely of the growing mode to second order in pertur¬ 
bations. This initial data is written in terms of a single 
function only, Smo- 


The expansion (85) has two conditions on its applica¬ 
bility: both the perturbative expansion and the derivative 
expansion must be valid. The expansion is constructed so 
as to carry equal weight to increasing the order of pertur¬ 
bations or the order of the derivative expansion. As such, 
this requires |5mo| 'C 1 and also \A5'^q/A + '^mol I'^mol 
(roughly speaking; of course, when Smo = O 7 this condition 
is difficult to satisfy, but violating it doesn’t necessarily 
mean that the derivative expansion fails to converge). 
As the solution (85) suppresses higher order terms more 
strongly as ^ — 00 , it is sufficient to ensure that these 

conditions are satisfied for the initial data at = 0. 


If the expansion fails to satisfy either of these criteria, 
this does not mean that the initial data is invalid; it may 
still satisfy the physical constraints above. All that this 
means is that the initial data no longer necessarily consists 
solely of the growing mode. It is not imperative that the 
decaying mode is nonexistent in the initial data; so long 
as the evolution starts with the perturbation sufficiently 
far outside the cosmological horizon, the decaying mode 
will have time to decay before the perturbation enters the 
horizon, which is when the typical black hole formation 
threshold is calculated. 


3. Linear Constraints 


We now write the above constraints in terms of the func¬ 
tion Smo, specialized to w = 1/3. Eor analytic purposes, 
we work to linear order only; it should be borne in mind 
that these constraints may be too strict when considering 
the effect of second order terms. Keep in mind that we 
are mostly interested in positive Smo, as this corresponds 
to an overdensity that can undergo gravitational collapse. 

Positivity of R requires 


(5 


/ 

mO 


< 


6(4 - Smo) 
A 


(93) 


while positivity of ifi requires <5^0 > — 1- 

To ensure that the various functions are even functions 
with zero derivative at the origin, we require <5^0 to also 
be an even function. Requiring R to be monotonically 
increasing yields the constraint 


8<5mO + < 


6(4 - Smo) 
A 


(94) 
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For m to be monotonically increasing {p > 0), we have 


0 < 1 + 6m0 + 


A 


^mO ■ 


(95) 


The requirement that F^ > 0 gives us the strongest 
constraint on the amplitude of 5mo- 

^ (96) 

Evidently, we cannot be even modestly nonlinear in Smo 
far enough outside the horizon of the initial data. In 
particular, for A > 2.5 (two and a half horizon scales from 
the origin in the initial data), fh is constrained to be less 
than 1.1, and so we are well within the linear regime using 
these variables. The only permitted nonlinear evolution 
that the initial data cannot capture is thus within the first 
couple of horizon spans. So long as the evolution begins 
with features of order the horizon scale or larger, this 
perturbative approach to initial data should be completely 
sufficient. 

For a black hole to be present in the initial data requires 
U to be negative, corresponding to 4 < 5mo- As this is 
well outside the linear regime, we do not expect any black 
holes to exist in the initial data constructed above. 

To match smoothly to FRW at large A simply requires 
i^mo —t 0 as A approaches the boundary of the computa¬ 
tional domain. Whenever the perturbative expansion ( |85[ ) 
is valid at the level at which is it truncated, the initial 
data consists of the growing mode only, by construction. 

The perturbative expansion requires |5mo| «C 1. It is 
harder to accurately state what the derivative expansion 
requires; a sufficient condition is |45(^g/A + |(5mol- 

One can also look to the derivative expansion terms in 
and require them to be less than 6^, which yields the 
constraint 




1 


A 




mO 


(97) 


To be entirely satisfied that the derivative expansion can 
be truncated with minimal error requires us to compare 
the next order term in the derivative expansion. Thank¬ 
fully, the derivative expansion terms are not difficult to 
compute. At first order in perturbations and up to second 
order in the derivative expansion, we have the following 
terms. 


p2« 


= 


Si = 


e 
30 

,35 

25^ 


S^S^Q e// 

T 


mO 


85; 


mO 


’mO 


A 


85; 


mO 


85) 


mO 


A2 


A3 


(98a) 

(98b) 

(98c) 


Assuming one can accurately calculate fourth order deriva¬ 
tives, it should be safe to truncate the expansion at second 
order in perturbation theory when the third order term 
here is small compared to the others. 


F. Relation to the Literature 


We now relate our results to other works in the litera¬ 
ture. 


1. Density and Mass Relationships 


We begin by writing the relationship between p 
in various ways. The exact relationship is given 
(50c), and can also be written as 


and fh 
by Eq. 


[{ARfm]' 

um 


(99) 


Using the result (62), this can be integrated to give 


m = I + ^^£[(Ai?)3]VA (100) 


where 5p = p — 1 is just the usual cosmological density 
contrast. The linearized version of these two relationships 
are 


_ [A35m]' _ A , nnil 

5p — A g 5^ (101) 

Sm^ ^ A^SpdA ( 102 ) 

with Sm = fn — 1. 

2. Fractional Mass Excess 


A number of papers compute the fractional mass ex¬ 
cess inside the cosmological horizon, given the (much- 
overworked) symbol 5. This quantity has been computed 
in different ways, with variations based on gauge and a 
choice on how the background mass is quantified. Carr’s 
original formulation PAIITS] used a uniform Hubble con¬ 
stant gauge with two shells of unequal density. A more 
recent approach by Harada et al. |56j revisited Carr’s 
approach. Furthermore, they conveniently represented 
their result in the gauge we have employed. 

The fractional mass excess associated with a perturba¬ 
tion is computed when the perturbation enters the horizon. 
In our language, the horizon is located at Ah = un¬ 
der background evolution; this is the comoving radius 
at which the comparison is performed The pertur¬ 
bation is said to enter the horizon when the density at 
the horizon is equal to the FRW density, followed by a 


® Note that the actual particle horizon will be slightly different due 
to nonlinear evolution. 
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sub-FRW underdensity to compensate for the overdensity 
inside the horizon. The time of horizon crossing is a 
gauge-dependent quantity, but within a given gauge, it is 
well-defined. 

The fractional mass excess S is then defined by 


0 = —1 
mb(AH) 

at this time, where mb{AH) is given by 

’mbiA) = —R^ijA^fjpb. 


(103) 


(104) 


In our language, the fractional mass excess is simply given 
by 


5 = Th{AH) — 1, (105) 

evaluated at the time of horizon crossing. 

Given an initial Smo, we can compute the expected 
fractional mass excess under the linear evolution de¬ 
scribed above. Let Ap^^rt be the comoving radius at which 
6p{Apert) = 0 in the initial data. The horizon will grow to 
this scale at time ^ = In(^pert)- Taking iirto account the 
growth of m during this time, the fractional mass excess 
will be approximated by 

<5 = SraoiApert)iApertf<'^-‘^^ , (106) 

although non-linearities typically become important be¬ 
fore the perturbation is completely inside the horizon. 

This result mirrors the result derived by Polnarev and 
Musco di], who defined 

S{A) = {m{A) - l)e2(“-b« (107) 


in our language, such that under linear evolution, <5 is 
constant. When evaluated at ^perti S = S, the fractional 
mass excess. 

Harada et al. [56] calculated that black holes are likely 
to form for S > Sc with 


5c 


5+ 3w \l + 3wj 


(108) 


This result is given in the gauge we employ here. For 
w = 1/3, this yields 6c « 0.4135. 

In order to extract the fractional mass excess from a 
numerical run, one should monitor the density p at the 
horizon. When it drops beneath one, read off 5 = to — I 
at that point in space and time. 

Using the fractional mass excess as an indicator of 
whether or not a black hole will form has some obvious 
deficiencies: if the density profile is such that the density 
contrast becomes negative briefly before becoming posi¬ 
tive again, this method will not capture the appropriate 
overdensity. Similarly, a long tail of a positive density 
contrast can also lead to erroneous results. In general. 


the dependence of black hole formation on the details of 
the density profile are not particularly well captured by 


the estimate in Eq. (108), although it works reasonably 
well as a rule of thumb. This is to be expected, as Harada 
et al. were more interested in correctly estimating the 
dependence of 6c on the equation of state than the initial 
profile. In order to overcome some of these deficiencies, 
Nakama et al. |49j have begun exploring more appropri¬ 
ate quantities for identifying black hole formation based 
on the density profile. 


3. Curvature Perturbation 


An approach to generating appropriate initial condi¬ 
tions based on curvature perturbations has been cham¬ 
pioned by Polnarev US] |37105] . This approach is based 
around the fact that 


F^ = 1 - KA^ (109) 

for a curved FRW universe. Accordingly, Polnarev and 
Musco 05] suggested letting 

= 1-K{A,t)A^ (110) 

and describing the dynamics of the curvature perturbation 
K (A, t). It was claimed that this picks out the growing 
mode; indeed, the growing mode is the particular solution 
of an ODE driven by the existence of the curvature per¬ 
turbation. The decaying mode exists as the homogeneous 
solution of the equation, which was neglected entirely. 
The method utilizes a time-dependent small parameter 
where 


e = 


2^.2 


a^r. 


(Ill) 


for some arbitrary, initially super-horizon comoving scale 
tq. It is useful to note that 




/o_„2(l-a)4 


( 112 ) 


Eventually, the various quantities 6\r are related to 
K{A,to) at first order in perturbations. The relation¬ 
ship between K and our variables comes from matching 
Eqs. (42el and (110). 

K{A, 0 = _ (j2'^ (^^3) 

Substituting to, U and R to linear order, we obtain 


K{A,0 = ^ (61^-261;) 


H 


Rjp 


(1 -b Of) 6^0 . 


(114) 

(115) 










16 


There are two ways to use K. The first is just to use it 
as a method to construct initial data. However, because 
of various approximations, when the actual curvature per¬ 
turbation is computed from initial data constructed using 
K, it doesn’t return exactly K. The second approach is 
to use K as an evolution variable to replace U. The issue 
with this approach is that K is related to instead of H, 
and as such is ambiguous whenever U may be negative. 

Despite these deficiencies, many aspects of this method 
were used as inspiration for our derivations, and we mostly 
arrive at the same results. Two major points that we 
changed were to construct completely dimensionless equa¬ 
tions, and to completely ignore the auxiliary variable K. 
Working directly with {7, m and R (or their tilded prox¬ 
ies) yields a more precise and clearer description of the 
system. In particular, our equations of motion are com¬ 
pletely absent of any scale (especially the scale ro which 
depends on the initial data), some auxiliary equations of 
motion have been eliminated completely, the growing and 
decaying modes have been distinctly identified, and all 
initial data is constructed based around the linear mass 
perturbation Smo- If desired, the curvature perturbation 
can always be computed. 


4- Bardeen Gauge Invariant Variables 

In order to connect to cosmological perturbation theory, 
it is useful to find the relation between our variables 
and the Bardeen gauge invariant variables [^. As the 
Bardeen variables are gauge invariant at linear order, 
here we work to linear order with our variables also. The 
relationships are derived in Appendix and reproduced 
here for convenienc43 

1 r°°-- 

$ = «'=/ ASmdA (117) 

2 JA 

As Sm grows as exp(2(l — a)^) in the linear regime, we 
see that <1’ and 'k are constant in this description, as we 
expect. We also find that the Bardeen variables are equal, 
as expected for a perfect fluid with no anisotropic shear 
stress. 

Heuristically, we can see that for super horizon scales, 
even very linear contributions to 6 m can cause nonlinear 
contributions to $ and 'k. This suggests that perturbation 
theory in 6 m and 6 u has a broader regime of validity than 
standard cosmological perturbation theory. We now show 
this more rigorously. 


The quantities and correspond to the Newtonian gauge 
metric 


ds^ = —(1 + 2^)dt^ + a^{l — 2'^)dx^ . (Hd) 


By using 6 p = 6 m + A 6 'm/i and integration by parts, 
we can write 


A 6 mdA = A^ 6 m ~ 3 / A5pdA 


(118) 


where the limit A^5m —t 0 as A —)■ oo by Eq. 
the squeeze theorem. Using this, we have 


(62) and 


$ = (^-A^ 6 m + 5j_ A 6 pdA^ . (119) 


Generally speaking, when 6 m is large and positive, the 
second term will be subdominant. To see this, consider 
some A beyond which 6 p < 0. Then 


A^6m = j I A\-6p)dA > 3 


A{- 6 p)dA (120) 


where the equality comes from integrating Eq. 
Thus, we can write 


( 101 ). 




1 

2 


A'^ 6 mO 


( 121 ) 


where we use the linear initial conditions. Requiring both 
$ > —1/2 and 'k < 1/2 to ensure a non-singular metric, 
we have 


6m0 < ^2 

. 1 
OmO > ^2 


( 122 ) 

(123) 


as conditions on the initial data (at linear order). These 
are very similar to the condition (96). In particular, this 
shows that the constraint (96) on 6 m 0 is not constraining 
the standard cosmological perturbation theory variables 
to be incredibly small at large radii. 


5. Connection to Inflation 

We now address the question of how to relate our 
evolution variables to a given profile generated by inflation. 
Our goal is to obtain an expression for 6 m 0 based on these 
initial conditions. 

The first step is to identify the initial time for our 
evolution. Consider that at the end of inflation, we know 
the Hubble scale Hi and the characteristic proper length 
scale of the perturbation, i?*. The ratio of the proper 
length scale of the perturbation to the horizon ratio, as a 
function of 5, is then 

n = (124) 

-^horizon vS ) 

where is the value of ^ at the end of inflation. We would 
like to start our evolution (^ = 0) when the perturbation 
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is still well outside the horizon, k ~ 5 — 10. This then 
sets to be 




1 

1 — a 


In 



(125) 


For typical large scale perturbations, will be somewhat 
large and negative. From Eq. (85), we see that both the 
perturbative and derivative expansion approximations 
become increasingly better for more negative As such, 
it should be sufficient to describe perturbations in this 
regime by using just the linear expansion. The goal of 
obtaining initial conditions is then successful once we 
obtain Smo- 


Let us assume that at the end of inflation, we have 
the gauge invariant metric perturbation (f> or T, and the 
gauge invariant velocity perturbation, V (see Appendix 
[§. From these, we can obtain Sm and 6u as 


6 


m 


Su 


A 

RrA 


(126) 

(127) 


where primes again refer to derivatives with respect to 
Knowing A provides no further information, as A = 5p, 
and Sm can be converted from one to the other. Note 
also that we cannot recover any information about Sr, 
which is gauge-dependent. This is unimportant; Sr has 
no effect on linear evolution, and rapidly converges to a 
solution that tracks Sm and Su appropriately. We should 
also point out that these formulas require ,V' —>■ 0 
sufficiently rapidly as A —>■ 0 (this is a requirement of 
spherical symmetry). 


Now that we know Sm and Sjj at the end of inflation, 
we need to write these in combinations that consist of 
the growing and decaying modes appropriately. From Eq. 
([72|), we have 


Sm = -b (128) 

Su = -b (129) 


where we make use of Eq. (69a). We are only interested 
in G, which we obtain through 


G= ^(5„-2<5c;)e2(“-i)«-. 
1 -b a 


(130) 


This provides a simple formulation to convert from 
Bardeen gauge invariant variables at the end of inflation 
to initial data for numerical evolution. 


^ Note that A is the coordinate system defined by ^4 = RrA, 
where Rr is the horizon scale at 5 = 0. One should be careful in 
translating a coordinate system at inflation into the coordinate 
system at the beginning of this evolution. 


V. HERNANDEZ-MISNER FORMALISM 


Though the Misner-Sharp formalism is very practical for 
the numerical study of spherically symmetric spacetimes, 
it breaks down if a singularity is formed. When simulating 
the formation of a primordial black hole, this makes it 
difficult to calculate its final mass. In order to deal with 
the formation of a singularity, Hernandez and Misner 
developed a formalism that uses radially outward traveling 
light rays to define a time coordinate m- The idea is 
to define a time coordinate based on the time at which 
an outgoing radial light ray starting from a given event 
reaches a distant observer. As time runs slower in a 
deep gravitational well, in this coordinate system, it takes 
infinite coordinate time for an event horizon to form, 
and so a singularity never forms. This formalism, which 
has typically been built on the work of Baumgarte et al. 
[36) . has been applied to primordial black hole formation 
in Refs. [3S1 0^ 031 03 03]. Our presentation of the 
formalism is somewhat simplified compared to previous 
implementations, as we eliminate the need for various 
thermodynamic quantities. 

Hernandez and Misner define a new coordinate u by 
the following differential relation. 

e'^du = e'^dt — e^^^dA ( 131 ) 


Here, is an integrating factor used to make du a perfect 
differential. Note that ii du = 0 then e'^dt = e^^'^dA, so 
along a path such that d6 = dcj) = 0 we have by Eq. ([^ 
that = 0. Eurthermore du = 0 implies that 


_ g0-A/2 

dt 


(132) 


Together, these relations show that u is an outgoing null 
coordinate. This coordinate system is known as “ob¬ 
server time coordinates”, as each u corresponds to what 
an observer at infinity sees looking towards the origin. 
Substituting Eq. (132) into the metric we find that 
the metric in this coordinate system is 


ds^ = -e^'^du^ - 2e^+^^^dudA + R^dn^ . (133) 


Given this metric, the equations of motion could be 
derived from Einstein’s equations. However, this is an 
incredibly tedious computation, and the same results can 
be arrived at by appropriately transforming quantities 
from the Misner-Sharp formalism. In order to carry this 
derivation out, we first define the following invariant op¬ 
erators. 


A = Dr = v^dfj, (134) 

Here, u^ = 0,0,0) is the comoving fluid four-velocity 

and = (0, 0, 0) in the Misner-Sharp coordinate 

system. By using the metric @ and the form of these 
operators in the Misner-Sharp coordinate system, we see 
that Dt is the derivative with respect to proper time at 
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a fixed spatial coordinate, and is the derivative with 
respect to proper distance at a fixed time and angle. 

We can write the Misner-Sharp equations in terms of 
these operators as follows. 


DtR = U 


(I35a) 

Dtm = -AttR^PU 


(I35b) 

^ fTDrP 

TO + AttR^PX 



P2 j 

(I35c) 

Drrn = AirpR? DrR 


(I35d) 

t~\2 -I tt2 2?Ti 

= I-\-U^ - 

R 


(I35e) 


Here we have used T = = D^R. We have omitted 

the equation for (j)' since this quantity does not appear in 
the Hernandez-Misner coordinate system. 

We now look at how various quantities transform when 
changing to the null coordinate. By looking at the trans¬ 
formation properties of the metric, we see that the cir¬ 
cumferential radius R is invariant under the coordinate 
transformation. Both P and p are scalar quantities that 
do not transform. The operators Dt and Dr are invariant 


by constr uction . Thus U is also invariant by Eq. (135a). 
Equation (135d I can be understood as a definitional equa¬ 
tion for to; it too is invariant under the coordinate trans¬ 
formation. When combined with the boundary condition 
m{A = 0) = 0, this defines m consistently across reference 
frames. 

Thus, to arrive at dynamical equations specific to the 
Hernandez-Misner coordinate system we can simply sub¬ 
stitute the appropriate expressions for Df and Dr in Eqs. 
(135). Through careful use of the chain rule, the deriva¬ 


tives become 


Here, overdots denote derivative with respect to u and 
primes denote derivatives with respect to A. In the density 
equation, we wrote P = p{w + Q), where we anticipate 
that Q will not depend on density; doing the same i n the 
U equation of motion leads to difficulties (see Section VH). 
Something to be cautious of in this coordinate system is 
that 


r = DrR = e-^l'^Rl - U (138) 

where the last term did not appear in the Misner-Sharp 
coordinates. Deriving the equation for is a little tricky; 
one must be careful with Drp. This is most easily calcu¬ 
lated by looping through the Misner-Sharp coordinates. 
Let 


and then use 


DrP = e-^l'^p - e~'^duP (139) 

e~'*duP = Dtp = e~’^dtp ■ (140) 


Using the result (20) followed by U'/R' = Drll/T to 
convert back to the null coordinate system yields the 
equation of motion for duU. 

The system of equations (137) is incomplete, as we 
require an equation for calculating ip. As dit is a perfect 
differential, we may write 


dtdAU = dAdtU. 


(141) 


and use the derivatives from Eq. (131) to obtain 


{Dt + Dr)^p = DrcP+]^DtX. (142) 


Dt = e-^du, Dr = e-^/^dA - e-^du ■ (136) 

Substituting these into the invariant form of the Misner- 
Sharp equations (135) and simplifying, we arrive at the 


Hernandez-Misner evolution equations. 
R = e^U 

rh = -c^AttR^PU 


(137a) 

(137b) 


[/ = -- 


p’/’ 


re-V2 


P' m + A-kRPP 


1 — w — Q 

J^e-’^Q + (u; + g) (2^+e 
P + P \ R 


P+P i?2 

ur 


p = 


e 


TO 


(137c) 

T + U 


47ri?2(r-(w-hQ)!/) 47ri?2i?'r - (w-K Q)t/ 

(137d) 


1 pt2 2 to 

= I -I- [72 —^ 


R 


e^/2 = 


R' 


T + U 


(I37e) 

(I37f) 


Next, note that we can rewrite Eq. |7b| as 

DrP 


Dr(j> = 


and Eq. (lib) as 


Dt\ = 


p + P 


2DrU 


(143) 


(144) 


The quantity Dk = Dt + Dr = e ^^"^Oa is a derivative 
along an outgoing null ray (as the derivative is taken with 
constant u). Combining all these results, we can write 
Eq. (142) as 


TO 


TDktp = DkU + ^ + 47r i?P 
R‘^ 


(145) 


where we employ Eq. (I35c). It turns out that this is not 


necessarily the best equation to evolve, however. Instead, 
apply Dk to (I37e) we obtain 


rDfeT = UDkU + 


D]^Rm 

P2 


Dkm 

R 


(146) 
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and substitute Dkin from (137dl, m/Er from (145), and 
DkR = [/ + r. After the dust settles, we obtain 

e'I’Dk (r + U)] = -AttR{p + P) (147) 

or alternatively, 

dAie-'^iT + U)] = -e^/^"’^47ri?(p + P). (148) 


When combined with an appropriate boundary condition, 
this closes the Hernandez-Misner system of equations. 

A number of boundary conditions from the Misner- 
Sharp formalism are invariant under the coordinate trans¬ 
formation. In particular, at the origin, i?(0) = 0, m(0) = 0 
and U{0) = 0. As previously, r(0) = 1. The conditions 
to match smoothly to FRW also remain the same as pre¬ 
viously: all local quantities must be equivalent to their 
FRW values, and the same mass matching condition must 
also be met. We discuss the outer boundary condition in 
the cosmological context below. 

Beware that R is not an odd function of A any more, as 
continuing the function to negative values of A continues 
to go back further in time. If instead the time coordinate 
increases again as A becomes negative, then R becomes 
even, but with a cusp at the origin. The implication of this 
is that quantities such as p which were even functions in 
the Misner-Sharp formalism are no longer even functions 
with zero derivative at the origin. 


The Hernandez-Misner formalism was originally for¬ 
mulated to describe a collapsing object surrounded by 
vacuum in an asymptotically flat spacetime. The time 
coordinate was taken to be the clock time of a static 
observer at future null infinity. As we do not have an 
asymptotically flat spacetime, such an observer doesn’t 
exist. Instead, we place an observer at the boundary of 
the computational domain. At this boundary, we know 
the quantities m, R, U, F, p and A, as well as the coor¬ 
dinates u and A. Let us also assume that we know the 
cosmic time coordinate that corresponds to this u and 
A (this can be found when initial data is constructed, 
see below). Then by the differential relation (131) with 
dA = 0, we have 


= gd-b _ 

dt 


(149) 


It is a convenient gauge choice to have u for this observer 
correspond to the observer’s cosmic time coordinate t 
also, which sets du/dt = 1 and thus tp = (j). Given the 
Misner-Sharp time, we can compute at the boundary 
through 


with 



(150) 


3q!^ 

S-kP 


(151) 
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Figure 2. Plot of null slices in Misner-Sharp coordinates for 
initial data from which a black hole forms. Slices of constant 
time in Misner-Sharp coordinates correspond to horizontal 
lines (constant ^). The way in which null slices avoid the 
formation of a black hole is evident, and corresponds to an 
essentially vanishing lapse. Plot data comes from the black 
hole evolution described in Section rVIIIl 


where we have used Eqs. (211, (351, 


and 


. When 


the density matches the appropriate FRW energy density, 

e'^ = 1. 


Let us now show that a horizon never forms in this 
coordinate system. The condition under which a trapped 
surface forms is 


U + r = DkR<0. (152) 


This relationship indicates that despite being evaluated 
along an outward-traveling null ray, the areal radius R of 
the ray is decreasing. From Eq. (1471 we see that 


Dk[e~^R>kRj = -e-^47rR(p + P) < 0 . 


(153) 


Thus, e~'^DkR decreases monotonically with increasing 
comoving radius. As e~'^Dj^R is positive for our observer 
at the boundary, it must then be positive everywhere. 
Thus, the condition for a trapped surface to form never 
eventuates, and a horizon never forms. We see how the 
null slicing condition avoids a black hole in Figure 
where we plot an example of null slices in Misner-Sharp 
coordinates. 


Musco et al. [IS] point out that as D]~R and can 
only reach zero together as u —>■ oo, care must be taken 
in the numerical evolution equations to ensure that this 
synchronization occurs, else negative values of D^R may 
result. The evolution equation (148) ensures the desired 
behavior regardless of numerical imprecision, while Eq. 
(1451 does not. 


The denominator F — r(;(l -|- Q)U in Eq. (137d) is of 
concern. This combination arises from 


DkiTi = 47rpi?^(F — (w + Q)U) 


Pb = 


(154) 
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where it is not guaranteed that F — (u> + Q)J7 is positive; 
at sufficiently large radii, the derivative of m along the 
outgoing null ray can become negative. However, as p is 
always positive, Di-m and T — (w + Q)U must vanish at 
the same event in spacetime. While a gridpoint is unlikely 
to exactly land on this point, numerical imprecision is 
likely to reduce accuracy in p in the vicinity of this region. 
For the background FRW evolution, the location of this 
crossover is i? = 1 /wH, at a radius a few times the horizon 
size (depending on w). So long as this radius is greater 
than the computational domain on the initial data for the 
Hernandez-Misner evolution, any issues related to this 
turnaround should be absent. 

For initial data generated as described below, we can 
compare the crossover radius to the size of the computa¬ 
tional domain for the FRW background. For the FRW 
background, the null ray’s position is given by 

1]. (155) 

When the null ray hits the boundary we have 

= 7W-W In ("1 + (1 - a) —^ . (156) 

(1-a) V a/ 

The crossover radius for the FRW background is 

(157) 

w 

which at becomes 



Demanding Ac S> A* yields the condition 

l»-t^^. (159) 

Thus, for all reasonable w, the crossover radius will be 
well outside the computational domain when the initial 
data is constructed, and this issue should not arise. We 
expect this result to still hold when using perturbed initial 
data, as the crossover radius is so much larger than the 
outer boundary, which we take to match to FRW (even if 
only asymptotically). 


A. Initial Data 


The biggest issue with this formalism is that it is un¬ 
common to posit initial data along a null slice. As such, it 
is typical to use the Misner-Sharp formalism to generate 
initial data, before converting to the Hernandez-Misner 
formalism. 

We suggest generating initial data for the Misner-Sharp 
formalism, and evolving this data forwards in time. In 
addition to evolving the physical variables, compute the 


trajectory of a null ray starting at the origin, and store 
the values of t/, m and R along the trajectory. If no 
black hole holes, then the Hernandez-Misner formalism is 
not required. If a black hole does form, then the domain 
can be truncated, excising the black hole, and evolution 
continued until the null ray reaches the outer boundary 
of the computational domain. 

As no inner boundary condition can be known for the 
truncated domain, we suggest repeatedly truncating the 
domain to remove unphysical artifacts. From experi¬ 
mentation, we found that setting a Dirichlet boundary 
condition of constant p at the inner boundary worked well 
to stabilize the boundary. Note that what occurs at the 
inner boundary cannot causally affect the evolution at 
the photon’s location, as the speed of sound is less than 
the speed of the photon. 

Once the null ray reaches the outer boundary, the data 
along this ray can be used as the initial conditions for the 
Hernandez-Misner formalism. This allows the evolution 
to be continued, and thus for the final black hole mass to 
be determined. 

In order to record the initial conditions for Hernandez- 
Misner evolution, a radial null geodesic needs to be 
tracked. This is described by du = 0, which is equiv¬ 
alent to 


dA 

dt 



(160) 


Writing this in terms of ^ and tilded variables, this is 


deA = ae^ - - . (161) 

d^AR) 

The initial condition is that the geodesic begins at the 
origin, A(^ = 0) = 0. 


VI. COSMOLOGICAL APPLICATION OF THE 
HERNANDEZ-MISNER FORMALISM 


In the Misner-Sharp formalism, we had much success 
in factoring the background cosmological evolution out of 
the evolution variables. This had the benefit of evolving 
dimensionless quantities that were all typically of 0(1), 
completely stable background evolution, and a simple 
manner in which to perform linear perturbation theory, 
which in turn led to an appropriate outer boundary con¬ 
dition. In the Hernandez-Misner formalism, it becomes 
more challenging to factor out the background evolution, 
as different values of A are located at different cosmic 
times t. In this section, we develop an approach that 
accomplishes the appropriate factorization. A similar ap¬ 
proach was suggested by Nakama et al. [IS], although the 
details of their implementation are somewhat opaque. 

We wish to evolve the variables U, m and R. After 
each timestep, the auxiliary variables F, p, P, and 
must also be computed. A further quantity describing 
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the cosmic time t{u, A) can also be evolved. This quantity 
obeys two relationships, based on Eq. (131). At fixed 
value of A, 


du 


while at fixed u, 


dt 

dA 


gA/2-,^ . 


(162) 


(163) 


As we did for the Misner-Sharp formalism, let us de¬ 


compose R, p, m and U as follows. 

R = RbR (164a) 

P = PbP (164b) 

m = -—pbR^rh (164c) 

U = HRU (164d) 

P = PbP = pbp{w + Q) (164e) 


All of the quantities appearing here are functions of u 
and A, where quantities that are known only explicitly as 
functions of t and A {Rb, Pb and H) are written as func¬ 
tions of t{u, A) and A. In practice, we use the previously 
defined ^ instead of t. We avoid ambiguous overdots from 
here on, instead expressing time derivatives explicitly. 

Let us assume that for a given m, we know R, U, rh 
and ^ as functions of A. As previously, we use A = RhA, 
and primes refer to derivatives with respect to A. We 
also define our time coordinate as u = Rhu in order to 
make it dimensionless. 

The equations of motion in these variables can be 
computed directly from the Hernandez-Misner equations. 
This is somewhat cumbersome, and so we relegate details 
to Appendix]^ At each timestep, a number of auxiliary 
variables need to be computed. The equations for these 
are as follows. 


p2 ^ g2(l-a){ _ jfi) 


(165a) 


P = 


= 


T + ARU 


AR m' — 2fh^' 


T-{w + Q)ARU 

3 aAk£,' + {Ak)' 


(165b) 

a^'Ak -k {Aky 

(165c) 

a^'{T + AkU) 


(165d) 


dA In (e-’^(r -k Ai?t/)) = -e' 


a — 1 


e^ARp 1 + w + Q 
f -k Am 1+w 


(165e) 


The last of these equations needs to be integrated in order 
to obtain e'^, using the boundary condition e'^ = e'^ at 
the outer boundary. 


Once the auxiliary variables are known, the evolution 
equations can be computed. 


duC = 

a 

duR = e'^-^R{U-e-‘^) 


dufh = ^ 


duU = -- 


e ‘^fh{l-\-w)— U{P+ 'm) 
V’-'f 


(166a) 

(166b) 

(166c) 


1 — w — Q 




a 


+ {w + Q)e^-^/‘^U' + ^ (ZQU 

AR{l + wPQ)\ 


w + Q 


AR{1 
- e-^d^Q) 


+ 3(1 -k w){U - e"'^) -k 

P 


(166d) 


Methods to evolve the U equation are discussed in Section 
IVIII below. 

In the Misner-Sharp formalism, we had the bound¬ 
ary condition = 0 at the origin for all tilded vari¬ 
ables. This translates to the coordinate invariant con¬ 
dition DrX = 0 at the origin, and is simply a result of 
the imposition of spherical symmetry. In the Hernandez- 
Misner formalism, this condition yields 

d^X = e^-^/^dAX (167) 


at the origin. 

The outer boundary is a little more complicated. Be¬ 
cause nonlinearities become important in the black hole 
formation process, our analysis for the outer boundary 
conditions in the Misner-Sharp coordinates is typically no 
longer valid. In particular, in practice, we found that both 
characteristics were directed inwards, and so no wave-like 
boundary condition is possible. The most stable outer 
boundary condition we found was to impose the Dirichlet 
boundary condition p = const, where the constant was 
determined by the initial conditions. While this tended 
to inject some energy into the system, there is no effect 
on the black hole mass, so long as the outer boundary is 
causally disconnected from the collapsing region. 


In order to match results from Misner-Sharp and 
Hernandez-Misner codes however (on a profile that doesn’t 
lead to black hole formation), it is useful to have the same 
boundary condition in place. In order to accomplish this, 
we can simply transform the boundary condition (66) 
we found in the Misner-Sharp formalism. In particular, 
derivatives can be transformed using Eq. (D17) and 


dA\t = dA\u-^^'"-^ d^\A ( 168 ) 


where we make use of Eq. (1361. The outer boundary 
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condition then becomes 


ae^+^-^d^6u = (j - I (C - 

- cs{6'u - e^/^-^duSu) - \S^ (169) 

where primes refer to derivatives in the Hernandez-Misner 
coordinate system. This can be rearranged to obtain 
duSu as follows. 


duSu = 


1 

Q,e<^+«-b - c^(,x/2-->p 

Cs _ 1 c| 

4 2 H 


X — r f)' 

2 ^ ^ ^s'^U 


(S' - e^^^~'^d-S ] - -(5 

\'^m ^ 


(170) 


This can be simplihed by using Eq. (Dll) to obtain 


duSu = 


a(l - Csf) 


11(^1+ c.e'-^e') 


1 ci 


2 A 


1 


+ l ^--^]SL-TSm- 


(171) 


Note that the (linearized) speed of sound is still Cg = 
e^/^/-\/l2, where ^ is evaluated at the outer boundary. 
This boundary condition works just as well as it does in 
the Misner-Sharp formalism. 

Obtaining initial data in this formalism is very straight¬ 
forward. Starting with the Misner-Sharp evolution, trace 
an outgoing photon geodesic as before, and save the val¬ 
ues for TO, U, R and Once the null ray hits the outer 
boundary, the data can be interpolated to the appropri¬ 
ate gridpoints, the initial value of u can be set using 
u = and the evolution can continue in the new 

coordinate system. 


A. Mass Extraction 


The black hole formation condition in the Hernandez- 
Misner formalism is the same as in the Misner-Sharp 
formalism: 


^ > 1. (172) 


Of course, this condition is never met; when it comes 
close to being met, the evolution is suppressed by the 
vanishing of the lapse. We thus need a method by which 
to extract the black hole mass. The method we describe 
here is based on that suggested by Baumgarte et al. [5S] . 
and demonstrated on an evolution that started with a 
Gaussian initial profile (see Section VIII for full details). 

In Figure]^ we plot the quantity 2m/R as a function 
of areal radius. We see that as the black hole gets closer 


1 
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Figure 3. 2m/R as a function of areal radius. The black hole 
formation condition is 2m/R > 1. In the inset, we plot the 
behavior near the formation condition. The lines in teal are 
earlier lines, and the lines in purple are taken at later times. 
We see that the formation condition is almost met at late 
times. The colors are used simply to help show the evolution. 
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Figure 4. The lapse e’^ as a function of areal radius. As 
with the previous figure, lines in teal are early, while lines in 
purple are later. The inset contains the same plot, but on 
a log scale. We see that as the black hole forms, the lapse 
plummets, essentially freezing evolution. A horizontal dashed 
black line is drawn at the 10“® level, and the corresponding 
radius Ro for mass extraction is read off from the vertical 
dashed black line. 


to forming, the curve approaches but never reaches the 
formation condition. As evolution continues, the line 
continues to approach 2m/R = 1, but never quite gets 
there. In our simulations, we found that 2m/R often grew 
to larger than 0.995 before we terminated the evolution. 

Figure]^ shows a plot of the lapse as a function of areal 
radius. We see that as the black hole forms, the lapse 
essentially vanishes at a reasonably well-defined radius. 
The inset shows that the lapse rapidly approaches zero in 
a logarithmic sense. 
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Figure 5. Mass enclosed in areal radius f? as a function of 
R (made dimensionless by Rh)- As with the previous plots, 
earlier timesteps are in teal, while later timesteps are in purple. 
We see that when the black hole forms, it sucks in all the nearby 
matter, and m becomes constant. If you extract the black hole 
mass at R — 6 Rh or R = IORh, you will find essentially the 
same mass. Off the right of the plot, m starts increasing again. 
A vertical dashed black line is drawn at l.li?o (extracted from 
the previous figure), and the mass contained in the black hole 
is read off from the corresponding horizontal dashed black 
line. 


In order to extract the black hole mass, we first hnd 
the radius at which the lapse starts to vanish. Program¬ 
matically, we chose to wait until the lapse at the origin 
satisfies < 10“^°. When this occurs, we find the areal 
radius i?o at which the lapse = 10“®. As the lapse is 
basically falling vertically at this point, it doesn’t partic¬ 
ularly matter if you choose 10“®, or a couple of orders of 
magnitude on either side. 

This radius i?o is the radius at which we wish to extract 
the black hole mass. In Figure we plot the function 
m{R) as a function of R. At R = Rq, we can read off 
the mass at which the lapse starts to vanish. However, 
note that the black hole has basically consumed all of 
the nearby matter. If we increase the radius slightly, we 
won’t change the mass of the black hole at all, but we 
will make sure to include everything in the elbow in the 
mass function. For this reason, we suggest reading off the 
mass at R = 1.1 i?o • 

Once the black hole has formed, it very slowly accretes 
mass, as evidenced by the slowly rising flat line of the black 
hole mass in Figure As this process is somewhat slow, 
it can basically be ignored for the purposes of estimating 
the final black hole mass. So long as one is consistent in 
the mass extraction procedure, the amount of accretion 
after initial formation should just be a small systematic 
error. 

In terms of our variables, mass is written as 

m = ^e-^/^R^A^rhRH ■ (173) 


So, once the radius you wish to read off from is known, 
that will need to be translated into a value for from 
which m and R will need to be extracted. 

Note that the mass is written in terms of the initial 
horizon radius scale. In order to compare this to more 
useful mass scales, we consider the mass within the horizon 
in an FRW universe (taking R = Rnit) = Rhs^)- 

R^horizoniC) ~ (^74) 

It is typical to compare the mass of the black hole with 
the mass contained in the horizon at the time when 
the perturbation enters the horizon. This time will need 
to be extracted from the Misner-Sharp evolution. Taking 
this ratio then eliminates Rh- 

For the example presented here, we find Rq = 5.79Rh- 
This radius corresponds to a mass of 2.89Rh, while l.li?o 
corresponds to a mass of 2.95Rh- The lines corresponding 
to these values are shown in Figures]^ andFrom looking 
at the Misner-Sharp data, we found that the perturbation 
entered the horizon at ^ = 2.47. The ratio of the black 
hole mass to the mass enclosed in the horizon at the time 
of horizon crossing is then 


'^black hole 
'^horizon crossing 


0.50. 


(175) 


Taking the black hole mass ten timesteps later in the 
evolution increased the mass we measured from 2.95Rh 
to 2.96Rh, or approximately 0.3%. The fractional mass 
excess from this evolution was approximately 0.43, just 
a little above the critical level of 0.4135 estimated by 
Harada et al. m- 

If you do not wish to go to the effort of getting the 
Hernandez-Misner formalism to run, you can always sim¬ 
ply estimate the mass of the black hole from the Misner- 
Sharp formalism; just leave the evolution to run until 
it dies from an inability to evolve further (due to the 
formation of a singularity). Then, find the outermost 
trapped surface, and read off the mass from there. 

This method tends to underestimate the black hole mass 
somewhat, as a fair amount of matter will fall through 
the horizon soon after its formation. For the example 
we have provided here, the mass we extracted before the 
simulation died was m = 2.25Rh, substantially less than 
the mass in the Hernandez-Misner formalism. We expect 
that the deficit is strongly dependent upon the initial 
density profile. 

In Figure]^ we plot time slices in Misner-Sharp coor¬ 
dinates as well as time slices in Hernandez-Misner coor¬ 
dinates. This figure shows where the Misner-Sharp mass 
estimate is extracted from, as well as the Hernandez- 
Misner mass estimate. This figure also shows that the 
Hernandez-Misner coordinates do an impressive job of 
hugging the event horizon, and that the Misner-Sharp 
coordinates can be used to probe a little bit inside the 
horizon. 
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In particular, if we look at the original null slice (at u = uq 
such that ^(0,?2o) = 0 and Eq. (178) holds), this yields 


dl-a)€(A„ 


,uo) — 1 1 _^4 

a 


(181) 


Looking towards from Eq. (145), we obtain 

= rH + AaH{a — 1)i9aC + ^ AH^a^ (182) 


where we make use of the Friedmann equation and writing 
aH in terms of Using our knowledge of from above, 
this becomes simply 


dA'tp = aH. 


(183) 


Figure 6. Time slices in both Misner-Sharp (teal) coordinates 
and Hernandez-Misner (purple) coordinates, plotted against 
areal radius. The Misner-Sharp simulation was run until the 
simulation could not proceed any further due to singularity 
formation at the origin. A red mark is placed at the outermost 
trapped surface on the final time slice, from which the mass of 
the black hole can be estimated. A second red mark is placed 
on the position at which we extracted the black hole mass 
from the Hernandez-Misner formalism. The data is the same 
as was used to generate Figure above. 


B. FRW Evolution 


For code debugging purposes, it is useful to know ana¬ 
lytic forms for all the variables in the case of pure FRW 
evolution. In particular, R = fh = p = U = e‘^ = lin 
FRW. It is straightforward to check that this is preserved 
under the evolution equations ( 166| ). Furthermore, in 
FRW, = a (the scale factor), and f = All 

that remains to be obtained then is ^{A,u) and 


We begin with From Eq. (Dll), we have 




gHRm 


a a 

This integrates straightforwardly to 




1 ~ cr T ^ 

- A + C 

a 


(176) 

(177) 


for some constant C. Our coordinate system is con¬ 
structed in such a way that on the outer boundary, 
t{Amax,u) = u- This implies that 


a 

Therefore, we have 

g(i-o)«A.^) ^ + /“A . ( 179 ) 

a \a/ 

By evaluating this at A = 0 and A = A^ax and taking 
the difference of the two results, we obtain 


^A-a)i{Ar„a^,u) _ g(l-a){(0,u) 


1 — a 


A 


max • 


(180) 


Writing this in terms of we obtain 


(184) 


where the prime denotes a derivative with respect to A 
as usual. We can use Eq. (179) to integrate this, yielding 


V' = 


1 — a 


In A — An 


1 — a Va 




+ C (185) 


for some constant C. We know at the outer boundary, 
e'^ = = 1, which sets the constant. We then obtain 


e 


b _ 



1 — a A — A„ 


a 


(^r 


a/(l — a) 


(186) 


For the initial null ray, we know from Eq. 
Uo = , which yields 


(178) that 


gbfT.uo) 


a -b (1 — a )A 
“b (1 f^^A^nax 


aj (1 —a) 


(187) 


Specializing to w = 1/3, these results become simply 
RO ~ (1 “f A^ax) 


= 1-b 


A-A 


max 




and 


gb(-4.«o) 


1 +A 

1 “b Ajjidx 


(188) 


(189) 


C. Domain Size and Computational Time 


In the Misner-Sharp formalism, doubling the size of the 
domain (i.e., doubling Amax) while maintaining the same 
resolution roughly doubles the computational time. In 
the Hernandez-Misner formalism, this is not the case, as 
we shall now show. 

Consider initial data in the Hernandez-Misner formal¬ 
ism where we specialize to w = 1/3, a = 1/2 for simplicity. 


a 
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From Eq. (1801, we can relate the value of u (which cor¬ 
responds to the value of ^ at the outer boundary) to the 
value of ^ at the origin, which we denote ^(0). 


W = ^ (190) 

Let us consider the change in u as the value of ^ at the 
origin progresses from ^(0) = 0 to ^(0) = 


Am = 


2 L 


e^-l + 2i™,,(e«/2_i) 


(191) 


We see that the change in u increases linearly with A^ax- 
This implies that if A^ax is doubled, then the computa¬ 
tional time to make a timestep equivalent to a timestep in 
the original domain also doubles, roughly speaking. This 
is independent of any increase in computational time from 
increasing the number of gridpoints in the domain. 

This problem is magnified when a black hole is forming, 
as the lapse near the origin becomes increasingly small. 
We thus see that taking a very large domain is compu¬ 
tationally expensive in the Hernandez-Misner formalism, 
much more so than the Misner-Sharp formalism. This is 
why we suggest using an appropriate boundary condition 
instead of steadily increasing the size of the domain in 
the Misner-Sharp formalism. 


VII. CODING ISSUES 



A 



In this section, we detail various issues pertaining to 
coding the formalism presented here. 


A. Evolution Scheme 

For a numerical evolution scheme, we suggest a finite 
difference method. The reason for this is that a lot of 
the interesting action tends to occur within a reasonably 
small radius, where a lot of gridpoints will be necessary, 
especially for profiles that are close to the threshold of 
black hole formation. 

The Misner-Sharp formalism is reasonably well-behaved. 
However, in the Hernandez-Misner formalism, we found 
that strong shocks tended to form in m (see Figure [^. 
These shocks were present in coordinate space, but were 
much less pronounced in physical space. If computing in 
the Hernandez-Misner formalism, we strongly suggest im¬ 
plementing an adaptive mesh refinement (AMR) scheme 
in order to evolve the system. An AMR scheme is re¬ 
quired to resolve the details of the system appropriately 
when such shocks occur, as gridpoints in comoving radius 
become very squashed and then rarified in the vicinity of 
a black hole (see Figure]^. 

In both formalisms, we performed time-stepping using 
an off-the-shelf adaptive time stepping Runge-Kutta (4,5) 


Figure 7. m as a function of comoving radius A (above) and 
areal radius R/Rh (below) in the Hernandez-Misner formal¬ 
ism. Data is taken from towards the end of the black hole 
evolution described in Section [VTII] These plots show the type 
of behavior one can except to see arising in black hole evolu¬ 
tions in Hernandez-Misner coordinates, demonstrating why 
finite difference methods with adaptive mesh refinement are 
necessary. A comparison of the plots shows that the spike in 
coordinate radius is not a physical shock: at around A ~ 2.8, 
the areal radius R{A) suddenly gets very large very quickly 
(see Figure]^. Thus, moving from A of ~ 2.8 to ~ 3 actually 
covers a large physical distance, which smooths out this feature 
somewhat in physical space. 


method. Starting from an initial profile, we could typi¬ 
cally extract a black hole mass in ^ 15-20 seconds on a 
consumer laptop. 

For gridpoints, a logarithmic grid spacing structure has 
been suggested in the literature (e.g., [H] and citations 
therein). Putting a lot of gridpoints near the origin and 
fewer at the outer boundary tends to put gridpoints where 
they are necessary to capture the dynamics of black hole 
formation. 

In terms of the outer boundary, we suggested above that 
having A^ax be too large would lead to a very slow evolu¬ 
tion. On the other hand, having A^ax be too small may 
lead to the outer boundary being in causal contact with 
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Figure 8. Plot of areal radius R as a function of comoving 
radius ^ in a black hole evolution (Hernandez-Misner formal¬ 
ism). The smoother lines come from earlier in the evolution, 
and the elbow develops around the same time that the black 
hole forms. The relationship between the comoving radius and 
the areal radius becomes progressively more squished when 
a black hole forms. This means that without an adaptive 
mesh refinement scheme, there will be very few grldpoints in 
A for a large change in R, which leads to poor sampling and 
inaccurate derivatives in this region. 


the collapsing region, in which case an imperfect boundary 
condition (particularly in Hernandez-Misner evolution) 
can lead to erroneous results. While the optimal size of 
Amax depends on the width of the initial perturbation, 
we found that a reasonable rule of thumb was to let Amax 
be 20 X larger than the initial perturbation. 

In the Hernandez-Misner formalism, we found that we 
had positive feedback for high-frequency noise in p. In 
particular, when a black hole came close to forming, p de¬ 
veloped a sawtooth instability, zig-zagging at the Nyquist 
frequency. We determined that this was due to our sam¬ 
pling scheme, and the instability was being amplified 
due to the presence of p' in the evolution equations. In 
order to compensate for this, when we calculated p in 
the Hernandez-Misner formalism, we applied a haystack 
smoothing filter over five gridpoints to p, which resolved 
the issue. This stabilization mechanism is a form of arti¬ 
ficial viscosity, separate from the viscosity we included in 
the equations of motion. 


B. Timestep Limits 


The Courant-Friedrichs-Lewy (CFL) condition [5D] pro¬ 
vides a limit to the size of a timestep that can be taken in 
a numerical evolution. Here, we provide the limit in both 
the Misner-Sharp and Hernandez-Misner formalisms |36j . 

In Misner-Sharp coordinates, the CFL condition re¬ 


quires 


A/2 A 4 

where Cs is the speed of sound, AH is the grid spacing, 
and At is the amount of time being stepped forwards, 
and we have used the speed of sound of a perfect fluid. 
This reduces to the condition 


At 

— < a — , _ 
to 


Writing this in terms of we have 


Ae = In ( 1 + ^ 


< In f 1 -b 


w 


(193) 


(194a) 

(194b) 


In the Hernandez-Misner formalism, we can use the 
discrete form of Eq. (1311 to obtain 

e‘^At = e’^Au-be^/2AH. (195) 

Eliminating At and solving for Au, we then obtain 

1 


Au < 


—-1 e^/2-’^AH. (196) 

w 


The CFL condition is a necessary but insufficient con¬ 
dition on stability. Further conditions on timestepping 
in the Misner-Sharp and Hernandez-Misner formalisms 
have been used in the literature |22L I36j , including lim¬ 
iting density changes to Ap/p < 0.02, radius changes to 
Ai?/i? < 0.005, and changes in the quantity l — 2m{R)/R 
to less than 0.1. Sometimes, these conditions can be 
much more stringent than the CFL condition. For our 
purposes, the adaptive time-stepper in our integration 
library proved sufficient, and we only demanded the CFL 
condition be met. 


C. Artificial Viscosity 


In our equations of motion, we included artificial vis¬ 
cosity from the very beginning, although we never gave it 
a form. We finally come to address this issue. 

Artificial viscosity is an additional pressure that is 
added into the system in order to mitigate the formation 
of shocks. The original implementations of the Misner- 
Sharp and Hernandez-Misner formalism included artificial 
viscosity in their equations of motion, and subsequent in¬ 
vestigations of supernovae collapse also utilized it. Papers 
in the literature using these formalisms for primordial col¬ 
lapse have sometimes implemented it, but documentation 
has been somewhat scarce. 

When discussing shocks, there are two distinct effects. 
Generally speaking, a shock is a discontinuity in a curve. 
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but can also simply describe very steep-sloping data. The 
second effect is a shockwave, which describes a moving 
shockfront or discontinuity. When shockwaves form, the 
characteristics tend to get distorted in such a manner as 
to guarantee that an actual discontinuity forms, which 
then propagates. This leads to all manner of difficulties 
when evaluating spatial derivatives. 

From our investigations, artificial viscosity, as it applies 
to the system that we have been investigating, is useful 
for smoothing out shockwaves, but much less effective 
at smoothing static “shocks”, which can simply be ad¬ 
dressed through increased local resolution. We found that 
outgoing shockwaves tended to form in the Misner-Sharp 
formalism when an overdensity failed to collapse into a 
black hole, and artificial viscosity was useful for smooth¬ 
ing these over. However, we never saw the formation 
of a shockwave in the Hernandez-Misner formalism; as 
a black hole forms, the vanishing lapse tends to ensure 
that everything freezes. We did see shocks forming, but 
artificial viscosity did nothing to assist with these. 

When shocks form in the Hernandez-Misner formalism, 
they tend to be a coordinate effect (see Figure for 
example). In our experience, these shocks can and should 
be suitably described through the implementation of an 
AMR scheme. 

Thus, our recommendation is to implement artificial 
viscosity in the Misner-Sharp formalism if you desire to 
describe outgoing shockwaves. We recommend ignoring 
artificial viscosity in the Hernandez-Misner formalism, 
noting that a few other authors have done likewise |44j[49]. 


1. Misner-Sharp Formalism 


The origin of the form of dissipation that we describe 
goes back to Neumann and Richtmyer |52) . which has 
been used in many evolutions of this kind [laiMiiiaETi- 
[5U] . although documentation is somewhat sparse. The 
basic idea is to identify where sharp discontinuities may 
occur, and artificially increase the pressure there in order 
to smear out the discontinuity over a few gridpoints. 

We investigated two triggers for identifying when to 
turn on artificial viscosity. The first condition, used by 
May and White [SB], is to trigger on dtp > 0. The 
second condition, used by Baumgarte et al. [35] (BST), 
instead triggers on BaU < 0. Both of these triggers 
seem reasonable: in the first instance, if the density is 
increasing, then something is collapsing. In the second 
instance, the idea is to identify shockwaves (moving in 
either direction) as requiring a negative velocity gradient. 

The suggested form of the viscosity term is 

_ J kp(AA 5^C/)2 when triggered . . 

\ 0 otherwise 



Figure 9. Density profiles from two Misner-Sharp evolutions 
from the same initial data. These evolutions did not form a 
black hole, but sent out a spherical wave of matter, which will 
tend to form a shockwave. In the purple evolution, artificial 
viscosity was turned on with k = 2, while in the yellow evolu¬ 
tion, artificial viscosity was turned off. The manner in which 
artificial viscosity has smoothed the shock is evident. 


take derivatives at constant t. 

In our description, we can construct Q as 


Q = K(AA)2e2(“-i)« 


Ba (aRu'^ 


The two triggering conditions become 


(198) 


A T^TJ' 

{AR)'U <--— (May and White) 

O 

{AR)'U < -ARU' (BST) 


(199a) 

(199b) 


where we use Eq. (20) for the first condition, and note 
that (AR)' > 0 by Eq. (89). Despite rather different 
origins, these conditions are very similar. Indeed, the two 
conditions are equivalent to taking the divergence of U in 
spherical polar coordinates {Ba{R^U) < 0) or in a linear 
sense {BaU < 0). When the velocities are converging 
instead of diverging, the viscosity triggers. 


We found that the BST trigger tended to work better 
than the May and White trigger, and we used it to smooth 
outgoing shockwaves in our Misner-Sharp evolutions. In 
particular, the smoothing provided by the artificial vis¬ 
cosity greatly increased the performance of our boundary 
condition. The effects of the artificial viscosity are evident 
in Figure]^ where we show data from evolutions with 
the same initial data, but with k = 2 and a = 0. In 
general, we found that k ~ 2 is a reasonable value of the 
coefficient. 


2. Hernandez-Misner Formalism 


where k is a dimensionless constant that controls the In the Hernandez-Misner formalism, we found that 
strength of the viscosity, AA is the grid spacing, and we artificial viscosity did not help much at all. The times 
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Figure 10. This plot shows three curves from late in the 
evolution of a black hole formation in the Hernandez-Misner 
formalism, plotted as a function of comoving radius.. The 
purple curve is U, which can be seen to have a negative 
derivative everywhere to the left of the dashed vertical line {U 
having a negative gradient is the trigger for artificial viscosity). 
The teal curve is m, which is the most shock-prone variable 
we evolve. We see that the shock is occurring outside the 
region in which artificial viscosity kicks in (to the left of the 
dashed vertical line), and so it does little to help smooth 
the shock. Furthermore, the lapse has essentially vanished 
to the left of the dashed curve, and so very little evolution 
is occurring anyway. The yellow curve represents R, which 
increases sharply around the point of the shock, causing a 
small distance in coordinate radius to correspond to a large 
distance in areal radius, which is the real cause of the shock. 
This leads to a dilution of gridpoints in physical space in 
this region, which can be addressed through the use of AMR. 
These curves have all been rescaled, so the vertical axis is 
unimportant. The horizontal black line shows zero on all three 
curves. 


when we thought that artificial viscosity may be helpful 
were when stationary shocks formed in the evolution 
variables. However, we discovered that artificial viscosity 
wasn’t triggering in the vicinity of these shocks. An 
example demonstrating this issue is described in Figure 


10 As the shocks that we saw forming were best addressed 


through the use of adaptive mesh refinement, we decided 
against the use of artificial viscosity in this formalism 
(setting Q = 0). However, for completeness’ sake, we now 
describe the implementation of artificial viscosity in this 
formalism. 


It turns out that transforming the artificial viscosity of 
the Misner-Sharp formalism into the null slices makes for 
a circular headache, due to Q depending on spatial deriva¬ 
tives in Misner-Sharp that correspond to both spatial and 
temporal derivatives in Hernandez-Misner. Instead, it is 
much easier to construct new conditions in the Hernandez- 
Misner coordinates. In particular, Baumgarte et al. [36] 
suggest triggering on 


SaU < 0 


( 200 ) 


in Hernandez-Misner coordinates, which for us is equiva¬ 
lent to 


Ai?Lf) < 0 . (201) 

Baumgarte et al.’s suggested form of the artificial viscosity 
is the same as in the Misner-Sharp formalism, but taking 
the derivative of U in the Hernandez-Misner coordinates. 
In our variables, this becomes 




( 202 ) 


This can be computed at every timestep. We again suggest 
K ~ 2 as a coefficient. 


Art ificial viscos ity appears in Eqs. (165bI, (165e), 
(166c I and (166d I in the Hernandez-Misner formalism. 
Two derivatives of artificial viscosity are required: OaQ, 
and duQ- The first of these is straightforward to evaluate. 
However, Q tends to trigger at a single gridpoint initially, 
which leads to somewhat unpleasant spatial derivatives. 
We suggest smoothing Q over a few gridpoints in order 
to assist with this derivative. 


The temporal derivative is somewhat harder. Baum¬ 
garte et al. suggest taking a one-sided first-order dif¬ 
ference (Q" — Q^~^)/Au to estimate this quantity. In 
most instances, this will be zero. The code will only be 
reduced to first order accuracy when shocks form and the 
triggering condition is met. 


VIII. IMPLEMENTATION EXAMPLES 


In this section, we discuss animations from simulations 
of the formalism that we have described in this papei|^ 
We have included still frames from the animations in 
this paper, but we believe that the animations provide a 
much better picture of the evolution than these frames 
can describe. For the benefit of those only interested in 
understanding the animations, we include a brief summary 
of the formalisms. 


A. Density Animations 


We start our animations with a 2D density map that 
evolves in time to form a black hole (“Black Hole Forma¬ 
tion” ). This animation requires no understanding of our 
formalism to interpret, and gives a good idea of how a 
black hole can form in the early universe. The expansion 
of the universe and the corresponding dilution of energy 


These animations may be found in a playlist on 
youtube at http://www. youtube. coin/playlist?list= 
PLqRPTTDr7yKTD9HExr7Fgk_RVrg5RwTau 
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Figure 11. Still frame taken from an animation describing the 
formation of a black hole in spherical symmetry. 


is evident as the map evolves, before the density collapses 
to form a black hole. The animation runs through twice; 
the first cycle shows the density using a linear scale, while 
the second cycle uses a logarithmic scale. A still frame 
from the end of the logarithmic version is shown in Figure 

El 


B. Misner-Sharp Formalism 

The rest of our animations are animated plots of our 
evolution variables. In order to understand them, it helps 
to know a little about the formalism we have used. 

We simulate the early universe using a constant equa¬ 
tion of state perfect fluid. For the purpose of these simula¬ 
tions, we take w = 1/3, corresponding to a radiation fluid. 
We use A as a comoving radius, and define a dimensionless 
comoving radius A = A/Rh, where Rh is the horizon 
radius at the time of the initial data. The Misner-Sharp 
formalism uses a cosmic time variable defined by appro¬ 
priate gauge conditions, and we evolve forwards using 
dimensionless time ^ defined by t = 

The Misner-Sharp formalism evolves three quantities 
forwards in time: the areal radius R{A,t), which is a 
gauge variable, the fluid coordinate velocity U{A,t), and 
the mass 

j-A 

miA,t)=4.TT pR^dARdA (203) 

Jo 

which can be thought of as the mass enclosed in the 
comoving radius A. The density p can be computed 
from m and R. For a number of reasons, it makes more 
sense to evolve quantities normalized by the values of the 
background cosmology, and so we define dimensionless 


tilded quantities i?, U, m and p as 


R = aAR = RbR 

(204a) 

U = HRU 

(204b) 

47r o 

TO = —pbRm 

O 

(204c) 

P = PbP 

(204d) 


where pb is the background (asymptotic) energy density. 
When evolving FRW, all tilded quantities are unity. It 
turns out that the initial background density po scales 
out of all the equations, and that all that is needed to 
specify initial data that consists only of a growing mode 
is an initial function for m = mg, from which initial data 
for Rq and Ug can be computed. 

For the next six animations, we demonstrate the evolu¬ 
tion of two Gaussian initial profiles for mg, taken a little 
over and under the threshold for black hole formation. In 
particular, we take 

rhg = 1 - 1 - Ke~"^ ( 205 ) 

with a = 2 Rh and ki = 0.173, K 2 = 0.17fQ We take our 
domain to be A^ax = 20Rh- In the animations, the teal 
curve represents Ki, which doesn’t form a black hole, and 
the purple curve represents K 2 , which does. 

After a black hole has formed, the singularity at its 
center rapidly forms, we can no longer evolve the simu¬ 
lation, and the purple curve in the animations vanishes. 
All of the plots displayed here are still frames from the 
animations, taken shortly after black hole formation. 

The animation “Density Profile” plots the density p/pg 
as a function of comoving radius A. The initial part of the 
animation looks incredibly flat, and shows the dilution of 
energy density due to the expansion of the universe. After 
some time, mass begins to accumulate near the origin and 
the black hole formation process begins. The density near 
the forming black hole remains many times higher than 
the background density, which continues to decay. As the 
purple curve forms a black hole, the density at the origin 
grows without bound. Our numerical simulation breaks 
down just before the singularity forms, and the purple 
curve vanishes at this point. When the teal curve fails 
to form a black hole, the extra energy is blasted out in a 
spherical wave. Note that the asymptotic (background) 
energy density drops at a constant rate, which occurs 
because we use a logarithmic time coordinate. The density 
profiles just after black hole formation are shown in Figure 

El 

The next animation “Density Perturbation (Comoving 
Radius)” shows the density profile relative to the back¬ 
ground density, which scales out the dilution due to the 


® Our first animation is of the energy density from the evolution 
using K 2 - 




30 



Figure 12. Density profile p/po as a function of comoving 
radius just after black hole formation. 


expansion of the universe. We see that the purple curve 
continues to increase in relative density until the black 
hole actually forms, while the teal curve reaches a peak of 
around a thousand times the background density before 
coming crashing down and releasing an outwards-moving 
density wave. The pressure when the teal curve peaks 
pushes the material out so rapidly that there is a strong 
rarefaction at the origin. The outgoing wave eventually 
forms a shockwave, which travels much faster than the 
speed of sound in the fluid. A still frame is shown in 
Figure 

The fourth animation “Density Perturbation (Physical 
Radius)” again shows the density profile relative to the 
background density, but this time plots it as a function 
of physical radius normalized by the horizon radius at 
the time (computed based on the background cosmology). 
Thus, R/Rh = 1 shows what is happening at the horizon 
scale at that time. Since the horizon grows with time, the 
horizontal scale corresponds to larger and larger physical 
distances as time progresses. We can clearly see when 
the perturbation has entered the horizon (when p drops 
below unity at R/Rh = !)■ Note that at the very end, the 
horizon scale grows larger than our domain of evolution, 
and so the outermost point we have data for recedes 
towards the left. It is interesting to compare the two 
plots in Figure by looking at features in these plots, 
we can see the difference between comoving radius and 
physical radius. In particular, the inner radius has been 
squashed somewhat, while the region where the teal curve 
lies above the purple curve has been extended. 

The fifth video “Velocity Perturbation” shows the nor¬ 
malized velocity perturbation as a function of comoving 
radius. Normalized velocities less than one are expand¬ 
ing slower than the Hubble rate, while negative veloci¬ 
ties are moving in the opposite direction of the Hubble 
flow, attracted into the growing overdensity. Initially, 
the velocity is essentially unity everywhere, indicating 
that matter is moving with the Hubble flow. We can 




R/Rnit) 

Figure 13. Density profile p as a function of comoving radius 
(above) and areal radius (below) just after black hole formation. 
The below curve has its x-axis scaled such that the horizon 
scale is at 1. 


clearly see when matter begins to fall towards the central 
overdensity. Eventually, the pressure of the overdensity 
dominates over the gravitational pull for the teal curve, 
and the velocity rises very sharply, explosively sending out 
a wave of matter (or photons, as the case may be). For the 
purple curve, the velocity continues to grow negative as 
matter continues to fall into the overdensity, until a black 
hole is formed. It is interesting to note a small apparent 
outflow {U > 1) from the black hole-forming region. This 
turns out to be a coordinate velocity effect; as the rela¬ 
tionship between physical velocity and coordinate velocity 
has become rather stretched in this vicinity, matter is 
still actually flowing inwards. The normalized velocity 
for the teal curve returns very rapidly to unity, despite 
being largely uninfluenced by the background cosmology. 
We believe that this is at least partly due to our gauge 
condition. A still frame can be seen in Figure [T^ 

Animations 6 and 7 (“Black Hole Detection (Comoving 
Radius)” and “Black Hole Detection (Physical Radius)”) 
demonstrate the evolution of the quantity 2m/R as a 
function of radius. The black hole formation condition 
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Figure 14. Velocity perturbation as a function of comoving 
radius just after black hole formation. 


is 2m/R = 1 (with G = c = 1; Eq. which signifies 

the formation of an apparent horizon (which must be 
contained within an event horizon). We see that after the 
purple curve crosses the yellow line, a singularity rapidly 
forms, and we cannot evolve it any further (the purple 
curve vanishes). There is a second condition for a black 
hole to form, namely that U must be negative at that 
point, which is why the curves can cross the yellow line at 
large radii without issue (the velocity there is roughly the 
Hubble flow, or U ^1). Still frames are shown in Figure 
|15[ Again, the difference between comoving coordinates 
and physical radius can be seen in the comparison between 
these two plots. 


C. Hernandez-Misner Formalism 



R/Rnit) 


Figure 15. Above: Plot of the black hole detection condition 
as a function of comoving radius. Below: The same plot, but 
as a function of normalized areal radius (such that 1 is the 
horizon scale at that time). A horizon has formed whe n a 
curve lies above the yellow line and U < 0 (cf. Figure 141. 


The Misner-Sharp formalism cannot evolve much past 
the formation of a black hole, due to the formation of 
a singularity. The Hernandez-Misner formalism aims to 
overcome this shortcoming by using a null time coordinate 
u. Data on a slice of constant u represents what an 
observer at the outer boundary of the computational 
domain sees at that instant, retarded by the distance 
that the light needs to travel to get to them. As no 
null rays ever leave an event horizon, a black hole never 
forms from the perspective of this observer. In particular, 
the lapse near the formation of a black hole approaches 
zero, which indicates that evolution is essentially frozen. 
Even though a black hole never forms, one can extract 
information about the black hole by looking at where the 
lapse becomes nonzero again. 

The Hernandez-Misner formalism evolves the same vari¬ 
ables as the Misner-Sharp formalism, just using a different 
coordinate system. The initial data for these variables has 
to come from a Misner-Sharp evolution, as constructing 
initial data on a null slice is somewhat unusual. A strong 


test of one’s code is to check to see that the evolution 
of a system is the same in the two different coordinate 
systems. This comparison is demonstrated in the eighth 
animation, “Two Coordinate Systems”. 

In this animation, we show a three-dimensional plot 
with comoving radius A on the cc-axis, logarithmic time ^ 
on the y-axis, and density perturbation p on the z-axis 
for an evolution that did not form a black hole. The red 
lines show timesteps from the Misner-Sharp evolution, 
where each time step is at constant The blue lines 
show timesteps from the Hernandez-Misner evolution, 
where each timestep is a null ray in the spacetime. The 
excellent agreement of the two surfaces shows that the 
evolution is in agreement between the two coordinate 
systems, including at the boundaries. A still frame from 
this animation is shown in Figure [T6| 

The rest of our animations are of the same purple black- 
hole-forming data that we showed in the Misner-Sharp 
animations above. 

Our next animation, “Density Profile (Null Slicing)”, is 
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Density Fluctuation Evolution 


Misner-Sharp 

Hernandez-Misner 




Figure 16. A plot of the density perturbation p as a function of 
A and ^ in both the Misner-Sharp (red) and Misner-Hernandez 
(blue) coordinate systems. 


Figure 18. This plot shows a frame from an animation of U as 
a function of areal radius, late in the formation of the black 
hole. As U is always negative, matter is sharply falling into 
the black hole. 



Figure 17. A plot of the physical density p as a function of 
areal radius. The large bump to the left has essentially been 
frozen in place. Together, the energy contained in this frozen 
bump constitutes the mass of the black hole. 


of the density profile in the Hernandez-Misner coordinates, 
plotting the physical density p as a function of areal radius. 
We see that the density prohle forms a large bump near 
the origin, as shown in Figure |17[ However, this bump 
is very much frozen in time by the vanishing lapse. In 
the Misner-Sharp coordinates, the density grows without 
bound at the origin; here, because we do not pierce the 
event horizon, it merely reaches a large density and stays 
fixed. As the evolution progresses, we see the energy 
density outside the black hole continues to fall off with 
cosmic expansion, though some matter accretes onto the 
black hole and freezes. Note that due to the null slicing, 
points on the curve at larger radius correspond to later 
times, and thus more expansion has occurred, decreasing 
the energy density. 


The animation “Velocity Perturbation (Null Slicing)” 
of the normalized velocity perturbation U as a function 
of areal radius shows the formation of a huge negative 
spike when the black hole forms. Figure shows a frame 
from late in the animation. We see that all matter is 
falling sharply inwards (recall that C/ = 1 is the Hubble 
flow). The animation also exhibits the freezing of the 
lapse. Starting from the origin, a trail of “frozen” velocity 
profile creeps outwards. 

Our final animation, “Lapse in the Null Slicing”, is 
that of the lapse as a function of areal radius. We do 
not include a still frame plot here, as the animation is 
simply an animated version of Figure]^ The animation 
clearly shows the lapse rapidly vanishing, which freezes 
the evolution of all variables. 


IX. CONCLUSIONS 

We have presented a comprehensive formalism for in¬ 
vestigating primordial black hole formation in full general 
relativity. The basic idea is to construct an initial density 
profile and evolve it forwards in time using a standard 
time-slicing until either a black hole forms or the overden¬ 
sity disperses. If a black hole forms, a second formalism 
that avoids the formation of a singularity is used to con¬ 
tinue the evolution and compute the mass of the resulting 
black hole. 

The foundation of the approach is the Misner-Sharp 
formalism, which we have adapted to cosmological appli¬ 
cations. We paid special attention to the outer boundary, 
where we constructed a non-reflecting boundary condition. 
For extracting information about black holes that form, 
we turned to the Hernandez-Misner formalism. Here, we 
adapted the system of equations to use variables more 
suited to cosmological and numerical evolution. 
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The issue of initial conditions was explored in detail, 
where we identified an expansion in two small parameters, 
and explicitly constructed initial conditions consisting 
only of the linear cosmological growing mode to second 
order, based on a single function Smo- The regime of 
validity of the expansion was investigated, including both 
perturbative convergence and physical constraints. We 
related these initial conditions to a variety of approaches 
in the literature, including Carr’s fractional mass excess, 
Polnarev’s initial curvature perturbation, Bardeen gauge 
invariant variables, and physical perturbation profiles at 
the end of inflation. 

The approach that has been developed here builds upon 
numerous other approaches in the literature. Most im¬ 
portantly, various ideas from the literature have been 
combined to form a clear and comprehensive picture. 
Numerically, there are numerous benefits to the enhance¬ 
ments we have incorporated, including completely stable 
cosmological background evolution, increased computa¬ 
tional efficiency, and increased numerical accuracy. The¬ 
oretically, we have investigated initial conditions in un¬ 
precedented detail, constructed a rigorous connection to 
inflation, introduced appropriate cosmological boundary 
conditions, and clarified the role of the initial curvature 
perturbation within this formalism. 

In addition to describing the formalism, we have con¬ 
structed a number of animations to describe what happens 
in primordial black hole formation and near misses. We 
hope that these animations will help spur interest in 
the subject, and provide a visual interpretation of what 
actually happens in primordial black hole formation. 

We plan to use the formalism developed here to inves¬ 
tigate properties of initial profiles that give rise to black 
holes and to quantify the resulting black hole masses that 
are produced. In the long term, we hope to run Monte 
Carlo simulations from inflationary initial profiles to iden¬ 
tify the number density and mass spectrum of primordial 
black holes formed from different models of inflation. 
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Appendix A: Stress-Energy Tensor 


In this appendix, we derive the stress-energy tensor 
for a spherically symmetric fluid, based on a statistical 
description of the system. This reproduces standard 
derivations, and is included for completeness. 

The description of a particle in a generic space-time 


is given by its position and momentum {a;“(A),p“(A)}, 
where an afhne parameter A is employed. The momentum 
is related to the position through 


P 


Ck 


da;“ 

~dA 


dx°‘ 


(Al) 


where the second equivalence is only defined for massive 
particles of mass m. A particle’s momentum obeys the 
mass shell equation, 

9°‘^PaPp = -rn^ . (A2) 


Note that given the spatial components of a particle’s 
four-momentum, this constraint determines the four- 
momentum time component, and thus only three of the 
four components of momentum are independent. 

Analogously to the definition in Special Relativity, we 
definite the number density of particles in phase space Af 
to be the function which describes the number of particles 
per unit volume in both position and momentum space. 
When integrated over any space-like hypersurface E and 
set of momenta TZ, this yields the total number of particles 
in E and TZ. 

N = J J dfj,{pi3)p^J\f{x^,pi3) (A3) 

s n 


Here, p°‘dT,a{x^) and dia{pp) are the invariant volume 
measures of E and TZ respectively. The 3-dimensional 
volume element dYia{x^) takes the form. 


dY,a{x^) 


dx^ dx" dx°' , „ , 

efiiycra - dadbdc 

oa ob oc 


(A4) 


where e^^a-a is the antisymmetric tensor with eoi 23 = +1, 
and we have parameterized E by (a, b, c) so that x“ = 
x“(a, 6, c). Since E is a 3-dimensional hypersurface in a 
4-dimensional manifold, its volume element is a vector 
pointing along the direction normal to it. In order to 
construct a scalar from dTia{x^) we must take its inner 
product with another four-vector. The only other true 
four-vector available in a statistical picture is the four- 
momentum , which we use to construct the invariant 
volume. 

The invariant momentum volume element is 


dp{pp) = 25 {p^pp + ni^) 0(p°) 


d^p 

7 ^' 


(AS) 


The delta function guarantees the momenta obey the 
constraint X2 and the step-function 0 picks out the 
positive energy solution. The factor of two arises from 
the fact that momentum is squared in the delta function. 
Note that the infinitesimals dp refer to momenta with 
lowered index (this arises because Pa is the conjugate 
momentum to x“). 


Note that in Minkowski spacetime on a constant time 
slice, we have dYia{x^) = (l,0)d^x, = {E,p), and 
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d^{p^) = so A3 


reduces to 


N = J (fix j (fip J\f{t,x,p) (A6) 


n 


as expected, where we have integrated over the delta 
function in momentum space. A number density per 
volume in space can be obtained by simply not integrating 
over the spatial hypersurface. 

We can take different moments of the number density Af 
by multiplying it by and integrating over momentum 
space. In particular, we can define the stress-energy tensor 
by the integral 

= J dp{p,)p^pdM{x\p,). (A7) 


The volume element d^p/yfi—g on momentum space 
then becomes the following. 


= e~’^dpo dpx dpy dpz 

= e~'^dpop‘^ sin Q dp dQ d^ (All) 


The range of p is from 0 to oo, while 0 < 0 < tt and 
0 < $ < 27r. The complete integration measure on 
momentum space then becomes 

dpip'') = 2S(^-^ sine dpded<^> 

(A12) 


Ignoring normalization, this is the only rank two tensor 
which we can construct from the invariants Af, p^, and 
dp{Pv) that is linear in A/”. One can be readily check that 
this object agrees with the usual stress-energy tensor by 
transforming to a local Lorentz frame. 

Let us now construct the stress-energy tensor in spher¬ 
ical coordinates. In order to compute the stress-energy 
tensor, we need to integrate over all possible momenta 
at a given position in space a;® = {r,9,(j)). To facilitate 
this, we construct an orthonormal Cartesian basis at a:*, 
oriented such that the x direction aligns with f, the y 
direction aligns with 9, and the z direction aligns with 
(j). We can then describe momenta at that point in space 
through the components {p^,p^,p^) in the orthonormal 
basis. For the metric described in Eq. (1^, these compo¬ 
nents are related to the momenta p^ = (p^, p^, p^, p'^) by 


where we take m = 0. 

We now calculate the stress-energy tensor. Inserting the 
integration measure (A12| into Eq. (A7) and integrating 
over the delta function, we arrive at 

T°‘9 = jP^ sine dpded^n’^n^Af (A13) 


where we have defined 


rp U— I 


= v°'lv = 


cos 0 sin 0 cos 4) sin 0 sin 4> 


=A/2 


R 


Rsin9 


(A14) 


In order for the distribution to be spherically symmetric, 
we require it to be a function of 


P = 


9 F 

P =R 


Pji 

y 


p-^ = 


p 


Rsin9 


(A8a) 

(A8b) 

(A8c) 


Af{x°',Pa) = Af{t,r,p,e) (A15) 


only, where cos 0 is related to the rotationally invari¬ 
ant dot product x’'pi. Note that p = \Jg^^PiPj is also 
rotationally invariant. 


Evaluating the integrals in Eq. 
component by component, we find 


(A13l over 0 and 4> 


We now write the components of the momentum in the 
orthonormal basis in spherical polar coordinates defined 

by 

(jf ,p^,p^) = (pcos0,psin0cos<l>,psin0sin$). (A9) 

Here, 0 denotes the angle away from the radial (x) axis, 
while $ is the azimuthal angle about the radial direction, 
and p^ = {p^Y + {P^)"^ + {p^Y- 
We can now construct Pr, pe and in terms of p, 0 
and 4). Begin by noting that Px = P^, thanks to the 
orthonormal basis. We can then construct the following 
relationships. 

Pr = e^^^Px = e^^^pcose (AlOa) 

Pe = Rpy = i?psin0cos4> (AlOb) 

p^ = R sin 9px = Rp sin 0 sin 4> (AlOc) 


= 27re-2'^y/dC (Al6a) 

= 2TTe-^-^/^ J CfdC (A16b) 

= 27re-^ J Cf dC (A16c) 

= ^ Ji^-OfdC (A16d) 

/(I - (A16e) 

sm 9 J 

with all other components vanishing, where we use ^ = 
cos 0 and define 

= J P^ffidp. (A17) 
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Each of these integrals involves different moments of the 
function / with respect to the variable C- 

If we assume that the particles in the spacetime in 
question are in thermal equilibrium, then at a given point 
in space, particles have equal probability to be going in 
any direction, and so / cannot depend on Under this 
assumption, the off-diagonal component vanishes, and if 
we set 


= (e-'^, 0,0,0) 

(A18a) 

P = 27r J fdC 

(A18b) 


(A18c) 

then we can write as 


= Pg°‘P + u°‘u^{p + P) 

(A19) 


which is a perfect fluid in comoving coordinates, and the 
desired result. Note that this also yields that the equation 
of state of radiation is w = 1/3. 

One might ask how non-interacting radiation (photons 
and neutrinos) can yield a fluid description. The key 
assumption that facilitated this was that of thermal equi¬ 
librium. So long as there are other species present in the 
universe to scatter radiation efficiently enough to main¬ 
tain thermal equilibrium, then the radiation bath behaves 
like a fluid. When neutrinos decouple from the bath and 
begin to freely stream, then the thermal equilibrium as¬ 
sumption will no longer hold, and the stress-energy tensor 
will develop an off-diagonal T*'’ component. At this stage, 
the photon bath may still be described by a perfect fluid, 
but the evolution of the neutrino species will have to be 
described by the full Boltzmann equation. Such a descrip¬ 
tion is beyond the scope of the formalism described in 
this paper. 


Appendix B: Method of Characteristics 

Here we use the method of characteristics on the lin¬ 
earized Misner-Sharp equations to find inward and out¬ 
ward moving components of the solutions. We want to 
find the inward moving component so that we can create 
a nonreflecting boundary at the outer edge of our compu¬ 
tational domain. We find the characteristic curves along 
which the solutions propagate, but the equations do not 
appear to permit a direct solution. For this appendix, 
we restrict our analysis to the case where w = 1/3. We 
also define cf = e^/12. Using these, we write the linear 
evolution equation, Eq. ( [70| ) as 

dlSm-ld^Sm-l5m-cl(^^+6':,^=0. (Bl) 

To use the method of characteristics, we need a first 
order linear equation, so we create a vector with the 


derivatives of dm and write the matrix equation; 


where 


dif + M ■ f = k, 



/O 0 0 \ 

M = 0 0 -c2 

VO -1 0 y 

^ ~ ^ (/i + h) + ■ 


(B2) 

(B3) 

(B4) 

(B5) 


To study the solutions of this equation, it is helpful 
for us to find the eigenvalues and eigenvectors of M. In 
this basis, we will still have a coupled system, but this 
will allow us to write down three equations where each 
has derivatives of only one the scalar functions appearing. 
The eigenvalues and eigenvectors are 



(B6a) 

(B6b) 

(B6c) 


We introduce Ui such that / = uiVi-\-U 2 V 2 from 

which we can identify 

5m = ui (B7a) 

d^5m = Cs{u 2 - Us) (B7b) 

6 m = U 2 + U 3 . (B7c) 

Substituting / into the equation of motion, we find 

yy Vi {d^iui) + Xiu)) = fc - . (BS) 

i i 


Note that because Cs depends on d^f has terms coming 
from the time dependence of the eigenvectors. Extracting 
the coefficient of each eigenvector from this equation, we 
obtain 


d^Ui = Cs (U2 - Us) 

1 2cs 

d(^u2 - Csu '2 = ^ ('“2 + M 3 ) 

1 2c 

d(^Us + CgU'^ = - ^{U2+ Us) ■ 


(B9a) 

(B9b) 

(B9c) 
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With the equations in this form, we may apply the 
method of characteristics. This method takes a first order 
PDE and computes curves called characteristics along 
which the equation becomes an ODE. We introduce A to 
parameterize these curves. 

d\ii = 1, = 0, dxui = Cs{u 2 - Ms) (BlOa) 

8x^2 = 1, dxA2 = -C, dxU2 = ^ ^ (U2 + Ms) 

4cs A 

(BlOb) 

dx^3 = 1 , dxA^ = c, dxus = ^ (m 2 + Ms) 

4cs A 

(BlOc) 

Here, Ai, and Ui are only functions of the parameter A. 

We can straightforwardly integrate the equations for 
and Ai, eliminating A. 


II 

(Blla) 

^ 2 (^ 2 ) = ^ 2 ( 0 ) — 2 cs 

(Bllb) 

^3(^3) = ^ 3 ( 0 ) + 2cs 

(Bile) 


Ai ( 0 ) are arbitrary constants giving the value of A on the 
slice ^ = 0. Thus we can see that the second family of 
solutions correspond to inward traveling waves, so if we 
wish to set a nonreflecting boundary, it is U 2 for which 
we need a boundary condition. 


Unfortunately, the equations cannot be solved explicitly. 
The condition we wish to satisfy is m^ = 0 for A > A^ax 
on the initial slice, ^ = 0 (assuming FRW initial conditions 
outside Amax)- As time passes and M 3 tracks outwards. 
Mi = 0 does not remain true beyond the boundary because 
M 3 sources the evolution of the other variables in Eqs. 
(BIO). Since beyond the boundary the values of mi and 
M 2 are initially zero and sourced only by M 3 , all of the 
information about their evolution should be encoded in 
the history of M 3 at the boundary A = A^ax- However, 
we could not find a way to use this knowledge to our 
adva ntage. Instead, we take the characteristic equations 
(B91 and pursue a local boundary condition, to reasonable 
success. 


Appendix C: Connection to Cosmological 
Perturbation Theory 


In this appendix, we connect the formalism developed in 
this paper to standard cosmological perturbation theory. 


1. Standard Cosmological Perturbation Theory 


We consider the perturbed FRW metric 
ds^ = — (1 + 2 (j))dt‘^ + 2ViB dx^dt 


(I + 2C)gi 


2ViVjE 


dx^dx^ (Cl) 


where gij is the background flat, time-independent three- 
dimensional metric on space, and Vi are covariant deriva¬ 
tives associated with this metric (we raise and lower in¬ 
dices on Vi with gij). We include only scalar perturba¬ 
tions, as that is sufficient for our purposes. Consider the 
transformation 

x^ ^ x^^ +i^^{t,xA (C2) 

with 

^0 = T, f = L* = V*L. (C3) 

Under this transformation, 

Sgf,^ - C^g^xu (C4) 


where g^i, is the background metric, and the Lie derivative 


is given by 

= i^dxQtiu + d^^^gxu + d^C^g^x ■ (C5) 

Evaluating the Lie derivative yields the following, where 
dots represent time derivatives. 

C^gu = -2T (C 6 ) 

C^gu = VAa^L-T) (C7) 

£ 55 ,, = 2a^{THg,, + V,V,L) (C 8 ) 

The metric scalars in Eq. (CD then transform in the 
following manner. 

(C9) 

B B + T - a^L (CIO) 

C-^C-HT (CIl) 

E^E-L (CI2) 


The Bardeen gauge invariant potentials are the fol¬ 
lowing combinations of these potentials. 


^ = cj) + dt{B-a'^E) (CI3) 

^ = -C - H{B - a^E) (CI4) 


We should also consider how the stress-energy tensor 
transforms under a gauge transformation. Let us split 
the stress-energy tensor into a background and perturbed 
component as 


(CI5) 

where 

= {pb + Pb)u'^ub. + Pb6i) (CI 6 ) 

with m(( being the background value of u^. Let us write 
the perturbed stress-energy tensor as 

= (Sp + 6P)u'{)Ubv + (pb + Pb){Su^Ub„ + u^Su^) 

+ SP6i) + (CI7) 







37 


where we again work only with scalar perturbations. We 


and use 


as 

= [l-py^v] 

(CIS) 

Up = [-(I -b (p), ^i{a^v + B)] 

(C19) 

no, = n% = 0 

(C20) 

/ -. - di - 

w.= v*V/--^v^ n. 

(C2I) 


The components of the stress-energy tensor 

are then 

= — (pf, -b dp) 

(C22) 

n =-(Pb + n)v*u 

(C23) 

= (p6 + Pfc)V,(a\ + n) 

(C24) 

T) = {Pb + 5P)5] + . 

(C25) 

Under a gauge transformation, the stress-energy tensor 

transforms as 


dT'", ^ dT'", - 

(C26) 

where the Lie derivative is given by 


- dx^T' 

V ■ (C27) 

Evaluating the Lie derivatives for defined above yields 

the following. 


C^f% = -Tpb 

(C28) 

= [pb + Pb)V^L 

(C29) 

=-iPb + Pb)^jT 

(C30) 

C^f) = TPb5] 

(C31) 

The scalars in the stress-energy tensor then transform in 

the following manner. 


Sp^ Sp- Tpb 

(C32) 

SP^SP- TPb 

(C33) 

V ^ V + L 

(C34) 

n 

(C35) 

A variety of gauge invariant quantities can be constructed 

for the matter perturbations. We focus on 

two of them. 

The first is the comoving density contrast A, defined by 

PfeA = dp -b Pb{a^v + B). 

(C36) 


The second, which we shall refer to as the gauge invariant 
velocity perturbation (as we are unaware of any standard 
name), is given by 


V = v + E. 


(C37) 


2. Connection to Misner-Sharp Formalism 


We now cast the Misner-Sharp metric ([^ into the above 
language. Writing 

R = aAR = aA{l + 6fi) 

A/2 _ 9aR _ 1 + + ASaSr 

r 1 -b dr 

where 

- ^drr?j 

and expanding the metric to first order in d/j, dr and (j), 
we obtain the following. 

ds^ = -(1 -b 2(j>)dt‘^ + 0^(1 -b 25R){dA^ -b A^dn^) 

-b a^{2AdA5R - 25r)dA'^ (C41) 


(C38) 

(C39) 

(C40) 


This metric contains three scalar perturbations, as ex¬ 
pected in a comoving gauge. No vector or tensor compo¬ 
nents are present. 


Looking back at the metric in Eq. (Cl I, let us take the 
metric gij to be the metric in spherical polar coordinates 
with radius A. We evaluate the quantity 


ViVjEdx^dx^ = d\EdA^ -b AdAEdn^ . (C42) 


Letting (p, C and E become functions of A and t alone, 
we can identify the coefficients in Eqs. (Cl) and (C41) to 
obtain 


Sr = C + 


OaE 

A 


dr = AOaC . 


(C43) 

(C44) 


The quantity (p is the same in both descriptions, and we 
have B = 0. 


We now relate our quantities to the Bardeen potentials 


(C13). In order to do so, we need to compute E, which 
in turn requires C. We thus start with C. 


C = -J^tdA' (C45) 

The limits of integration are set by the requirement that 
all perturbations vanish at infinity. In order to compute 
C, we first compute 


deC = - 


^dA' 

A' 


(C46) 


where 


Sjdr = 2 (a - I)dr + (^d^6u - 


= -A^—OaSp 


(C47) 
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simplify this expression. We then have 


and we have used Eqs. ( [M| and a(l + 3w) = 2(1 — a) to 

(C48) 


d^C= - —5p. 


Converting from ^ to t, we obtain 




(C49) 


Next, we look towards E. Taking our expression for Sr 
and integrating as we did for C, we obtain 


E = Rjj A'{C - Sii)dA'. 


(C50) 


We t hen c omp ute d^E as a step towards E, using Eqs. 
(68a) and (68b) to simplify the result. 


d^E = -Rjj / aA'SudA' 


E = -HR]j / A'SudA' 
Ja 


The first gauge invariant variable is then 


T = / A'SmdA' 

2 Ja 


(C51) 

(C52) 

(C53) 


where we use the definition of (5r (C40). 


Eor the second gauge invariant variable, we need 
dt{aJE). As before, it is simplest to calculate d^{aJE) 
and then convert the result. 


d^{a^E) = ^a^HRl J_ A' S^dA' - ^Sp 

dt(a^E) = A' S^dA' - ^Sp 


(C54) 

(C55) 


Here, we’ve used Eq. (68e) to substitute d^Su- The 
second gauge invariant variable is then 

1 


$ = / A! 5mdA' = A! 

2 Ja 


(C56) 


where we substitute </> from Eq. ( |68a[ ). We could have 
anticipated this result; a perfect fluid has no anisotropic 
shear stress, and so the two Newtonian potentials, and 
thus the Bardeen gauge invariant variables, are the same. 
Note that as our linear solution for Sm grows as 
$ and T are invariant to first order in our expansion. 

The gauge invariant density contrast is simply A = dp, 
thanks to the comoving gauge and diagonal metric. The 
gauge invariant velocity perturbation is 


V = E = -RHe-^ / A'dudA'. 
Ja 


(C57) 


Appendix D: Derivation of Hernandez-Misner 
Results 


In this appendix, we include some derivations relating 
to the Hernandez-Misner results. 


1. Hernandez-Misner 

In the Hernandez-Misner formalism, the equation of 
motion for U needs to be taken carefully. We begin with 
Eq. (I35c). Substituting for the derivatives in Eq. (136), 


we obtain 
e-^d„U = - 


r , , TO-I-47ri?^P\ 

(D,P-AP) + -5 -,-- \ ■ 


p + P 




(Dl) 


We can write 

DtP = Dt[piw + Q)] = Dtpiw + Q)+ pDtQ. (D2) 

Inserting this back into the equation of motion, we obtain 
the following. 


-^duU = - 


T 


{DrP - Dtp{w + Q) - pDtQ) 


p + P 

TO -I- AttR^P 

A 


(D3) 


We then insert Dtp = e ‘^dtp using Eq. (20), being 


careful with the Misner-Sharp spatial derivatives of U 
and R. 

e-*ajj = 

p + P R‘‘ p + P 

-{w + Q) ( 2 ^ + e-^l'^U' - e-^duU^ 

(D4) 

Primes here refer to derivatives with respect to A in 
Hernandez-Misner coordinates. This can be solved for 
duU as follows. 


duU = -- 


oi> 


1 — w — Q 


re -^/2 


P' m + AttR^P 


p + P 


i ?2 


-^e-^d^Q + {w + Q)(2^+ 
p + P \ P 


(D5) 


2. Cosmological Hernandez-Misner 


We now derive the cosmological equations of motion 
for the Hernandez-Misner formalism. In particular, we 
use the definitions (164) in the equations of motion (137) 
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and (148). We assume that at any given time u = Rru, 


we know U, rh, R and ^ as functions of A. 

The first step is to compute F. As previously, it is 
useful to work with f = T/HaRn- 


p 2 ^ g 2 (l-a)« a^R^{U^ - m) 


Next, we construct p from Eq. (137d). 


P = 


m 


T + ARU 


AirptR'^R' T-{w + Q)ARU 

To expand this further, we need 

R' ,, 1 R' 

R-^^^A^h 

and 


m' = AnpbR^rh 


(D6) 


(D7) 


(D8) 


V 

,, 1 

R' 

1 fh 

a — - 


H —^ 

+ — —^ 

_V 

A 

R 

3 m 


(D9) 


The end result is 

f + Am 

P = 


T-{w + Q)ARU 


fa + 


AR fh! — 2fh^' 

aAk^' + {Aky 


(DIO) 


Knowing p, we can compute P = p{w + Q) (given a form 
for Q). 

We now turn our attention to From Eq. (163), 
we find 


„A /2 _ 




From Eq. (137f|, we have 


e -^/2 ^ g 


.r + AR{U - e-^) 

{A^' 


(Dll) 


(D12) 


As we do not yet know e"^, we can use these equations to 
eliminate the unknown and solve for e‘^, and then 
compute using Eq. (Dll). 


_ a^'AR + {Ary 
~ a^'{r +Am) 


(D13) 


The final auxiliary variable we need to compute is . 
From (1481, we have 


dA[e-'>’iT + ARU)] 

= 


e'^ARp ^^^'^^^ + (a - l)(r + ARU) 


1 + w 


(D14) 


subject to the boundary condition that = e‘^ at the 
outer boundary. To integrate this equation, let y = 
g-b(f _|_ ARU), and rewrite it as 


'9A(lnx) = -C' 


e'f’ARp 1 + w -\- Q 

f + Am 1 + w 


+ a — 1 


(D15) 


which can be integrated straightforwardly. 

The next step is to compute the evolution equations for 
R, U, fh and For fh and R, the results can be found 
by changing coordinates in the Misner-Sharp evolution 
equations in these coordinates, Eqs. (50). In particular. 


DtX = e-^dtX = e-*duX . 


(D16) 


Writing this in terms of derivatives with respect to u and 
X, we have 

d^X = ae^+^-^duX . (D17) 

Applying this to R and fh, we obtain the following. 

(D18) 

duR = e'<’-^R{U-e-‘^) (D19) 


dui = , 

a 


dufh = 3e’^ ^ 


e ‘^fh(\ + w) — U{P + fh) 


(D20) 


Note that when the lapse —)■ 0 near horizon formation, 
all quantities are essentially frozen; we find the same 
behavior in the U evolution equation too. 

To obtain the equation of motion for U requires a 
bit more work. Start with Eq.(50f|, transform the 9^ 
derivative as above, and write the spatial derivatives in 
terms of the invariant derivative operators. 


.U 


DrP 


e^-^duU = e-‘^ - 

a ARDr{AR){p + P) 

1 


- ^ (2(7^+m + 3P) (D21) 


Now, we can write 


Dr{AR) = ^ (AR) 


Rh 

= Ht 


(D22) 

(D23) 


where we make use of the result for duR above and Eq. 
(D12|. The other derivative we need to evaluate is 


DrP = (iC + Q)DrP + pDrQ ■ 
The first term here can be written 

DrP = - Dtp. 

Hh 


(D24) 

(D25) 


To evaluate the final term here, we turn to Eq. (20) in 
the form 


, ,, U DrU 

Ap--(p + P)( 2 - + — 


(D26) 
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Writing this in terms of tilded quantities, we have 

Ap = 377e-^l + w)p - H{p + P)(w+ 

\ Dr{AR)) 

(D27) 


This equation looks like it diverges as A —>■ 0. In order 
for it not to diverge, we require 


where we use the continuity equation in the form DtPb = 
—e~‘^3Hpf,{l + w). Combining all these results, we can 
write DrP as 


DrP = {w + Q) 


pDyQ 

Rh ^ w + Q 


RA 


+ (p + P)^{e-^/^U' - e-^d,,U) 


3Hp 


{1 + w){U - e-^) + QU 


(D28) 


where primes once again refer to derivatives in the 
Hernandez-Misner coordinate system. 


We can now combine these results in Eq. 


(D21I. 


lim 4 + 3(1 + w){U - + 3QU 

A-s-O [ p 


= 0 . 
(D31) 


This is a sufficient condition, as this quantity is just being 
divided by A. (Note that DjiQ = 0 at the origin by virtue 
of being an even function in Misner-Sharp coordinates.) 
In order to show this, we write p in terms of Dr and Dt 
once again. 


,t/ 1 

a 

f (w + Q) 
ARH{p + P) 


e^-'^duil = - + 3p) 

a 2 \ ) 

'e-V 2 


Rh 


pDrQ 

w + Q 


+ 3Hp 


[I + w){U - e-*) + QU 


p' = RHe^/^{Dt+Dr)p (D32) 

As p is an even function in Misner-Sharp coordinates, 
DrP vanishes at the origin, and we can write 


-{w + Q){e^-^/^U' - e^-'>’duU) 
We can finally solve this for duU. 


duU= - 




IP 


1 — w — Q 

+ {w + 

_ f{w + Q) 


3P] +U^-e-^- 
a 


(D29) 


U 




DrQ 


AR{\ -|- u; -|- Q) 

+ 3{l + w){U - e-’t’) + 3QU 


p H{w + Q) 


(D30) 


p\A = 0) = [e"'^(l -h w)p -{p + P)u] (D33) 


using Eq. (D27) from above. Inserting this in Eq. (D31|, 
it is straightforward to show that the limit vanishes. We 
thankfully do not need to evaluate d^U at the origin, as 
there is a boundary condition there, but it is good to 
know that the equation of motion is not divergent as one 
approaches the origin. 
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